High-performance linear-system solver for shielding current analysis in cracked HTS film

Abstract

A high-performance method is developed for solving a linear system in the shielding current analysis of a high-temperature superconducting (HTS) film containing cracks. When the shielding current density is calculated in the HTS film, a linear system of special type has to be solved at each iteration cycle of the Newton method. If GMRES is applied to the solution of the linear system, its convergence property will be remarkably degraded with increasing crack size. In order to resolve this problem, other variables than corrections of the current vector potential are all eliminated from the linear system. Consequently, the residual history of GMRES is hardly affected by a crack size so that its convergence property is improved remarkably.

Keywords

Finite element method high-temperature superconductor linear-system solver Newton method

1. Introduction

Since evaluation of the shielding current density is indispensable in designing engineering applications of a high-temperature superconducting (HTS) film, several numerical methods [1–4] have been proposed to anlyze the time evolution of the shielding current density. After discretized with respect to time, an initial-boundary-value problem of the shielding current density is transformed to a nonlinear boundary-value problem at each time step. However, the solution of the nonlinear problem by the Newton method is extremely time-consuming especially for the case with a cracked HTS film. This is mainly because the convergence property of the linear-system solver, GMRES, is remarkably degraded for this case.

The purpose of the present study is to develop a high-performance method for solving a linear system in the shielding current analysis of a cracked HTS film.

2. Shielding current analysis

A cracked HTS film of thickness b is exposed to the time-varying magnetic field. We first assume that it contains m cracks and that it has the same cross section Ω over the thickness. By taking its thickness direction as z-direction and choosing its centroid as the origin, we use the Cartesian coordinate system 〈O : e _x, e _y, e _z〉. In the following, j and E are a shielding current density and an electric field, respectively. Also, 𝜇₀ and B denote the permeability of vacuum and an applied magnetic flux density, respectively. Furthermore, n and t denote an outward normal unit vector and a tangent unit vector on the boundary $\partial$ Ω of Ω, respectively, and x and x ^′ are position vectors of two points in the xy plane. It must be noted here that the boundary $\partial$ Ω is composed of not only the outer boundary C ₀ but also crack surfaces, C ₁, C ₂, …, C _m.

Under the thin-plate approximation, there exists a scalar function T ( x , t) such that j = (2∕b)∇ × (T e _z), and its time evolution is governed by the following integro-differential equation [1,2,5]: $\begin{eqnarray}\displaystyle {\mu}_{0}\frac{\partial }{\partial t}({\hat{W}}T)+({\rm\nabla}\times \boldsymbol{E})\cdot \boldsymbol{e}_{z}=-\frac{\partial }{\partial t}\langle \boldsymbol{B}\cdot \boldsymbol{e}_{z}\rangle . & & \displaystyle\end{eqnarray}$ (1) Here, 〈〉 denotes an average operator over the thickness and ${\hat{W}}$ is an operator given by $\begin{eqnarray}{\hat{W}}T\equiv \frac{2T(\boldsymbol{x},t)}{b}+\iint _{{\rm\Omega}}Q(|\boldsymbol{x}-\boldsymbol{x}^{\prime }|)∼T(\boldsymbol{x}^{\prime },t)∼d^{2}\boldsymbol{x}^{\prime },\end{eqnarray}$ where Q (r) = −(πb ²)⁻¹[r ⁻¹ − (r ² + b ²)^−1∕2].

The initial and boundary conditions to (1) are assumed as follows: T = 0 at t = 0 and j ⋅ n = 0 on $\partial$ Ω. For the J–E constitutive relation, the following power law [3–8] is assumed: $\begin{eqnarray}\displaystyle \boldsymbol{E}=E_{\text{C}}\left(\frac{|\boldsymbol{j}|}{j_{\text{C}}}\right)^{N}\frac{\boldsymbol{j}}{|\boldsymbol{j}|}, & & \displaystyle\end{eqnarray}$ (2) where j _C and E _C denote the critical current density and the critical electric field, respectively, and N is a positive constant. By solving (1) with the initial and boundary conditions, we can analyze the time evolution of the shielding current density j .

The above analysis is applied to the scanning permanent magnet method (SPMM) [9] that is one of contactless methods for measuring spatial distributions of j _C in HTS films. In the SPMM, a cylindrical permanent magnet of radius R and height H is moved along the surface of an HTS film and, simultaneously, an electromagnetic force F _z acting on the film is monitored (see Fig. 1). During the movement of the magnet, the distance L between the magnet bottom and the film surface is kept constant. In the following, the cross section Ω of the film is assumed as a rectangle of width w and length Aw. Furthermore, the longitudinal direction of the film is taken as x-axis. Also, the symmetry axis of the magnet is denoted by (x, y) = (x _A, y _A), and its movement is assumed as x _A = vt − Aw∕2 and y _A = const. Here, v is a scanning speed. Moreover, the film is assumed to contain a single crack whose cross section is a line segment connecting two points, ( ± L _c∕2, 0 mm), in the xy plane. Throughout the present study, the physical and geometrical parameters are fixed as follows: j _C = 1 MA∕cm², E _C = 1 mV/m, N = 20, b = 1 μm, w = 12 mm, A = 11, m = 1, R = 0.8 mm, H = 2 mm, L = 0.5 mm, y _A = 0 mm, and v = 10 cm/s.

Fig. 1.

A schematic view of the scanning permanent magnet method.

3. Virtual-voltage method

The authors developed the virtual-voltage method [4–6] for accurately solving the initial-boundary-value problem of (1). In this section, the basic idea of the virtual-voltage method is explained briefly. In the following, a superscript (j) denotes a value at time t = jΔt, where Δt is the time step size. In addition, the voltage around C _i and its numerically evaluated value are denoted by h _i( E ) and N _i( E ), respectively. For the numerical integration of h _i( E ) ≡ ∮_{C
_i} E ⋅ t ds, the Gauss-Legendre quadrature is employed. Here, s is an arclength along $\partial$ Ω.

In the virtual-voltage method, T ^(j)( x ) is determined as a solution T of the following nonlinear boundary-value problem: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle G(T)\equiv {\mu}_{0}{\hat{W}}T+{\rm\Delta}t∼\boldsymbol{e}_{z}\cdot ({\rm\nabla}\times \boldsymbol{E})-u=0∼\text{in}∼{\rm\Omega},\end{eqnarray}$ (3) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle T\in H(\bar{{\rm\Omega}})\equiv \left\{w(\boldsymbol{x}):w=0∼\text{on}∼C_{0},∼\frac{\partial w}{\partial s}=0∼\text{on}∼C_{1},C_{2},\ldots ,C_{m}\right\},\end{eqnarray}$ (4) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle h_{i}(\boldsymbol{E})={\phi}_{i}\quad (i=1,2,\ldots ,m),\end{eqnarray}$ (5) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle N_{i}(\boldsymbol{E})=0\quad (i=1,2,\ldots ,m),\end{eqnarray}$ (6) where u ≡ 𝜇₀ ŴT ^(j−1) −〈( B ^(j) − B ^(j−1)) ⋅ e _z〉. Note that the above nonlinear boundary-value problem is solved for an unknown function T and unknown constants $\{{\phi}_{i}\}_{i=∼1}^{m}$ .

If the Newton method is applied to (3)–(6), we get the following linear boundary-value problem: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle {\delta}G=-G(T),\end{eqnarray}$ (7) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle {\delta}T\in H(\bar{{\rm\Omega}}),\end{eqnarray}$ (8) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle {\delta}h_{i}-{\delta}{\phi}_{i}=-[h_{i}(\boldsymbol{E})-{\phi}_{i}]\quad (i=1,2,\ldots ,m),\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle {\delta}N_{i}=-N_{i}(\boldsymbol{E})\quad (i=1,2,\ldots ,m).\end{eqnarray}$ (10) Here, δT and δ𝜙_i are corrections of T and 𝜙_i, respectively, whereas δG, δh _i and δN _i denote Fréchet derivatives of G (T), h _i( E ) and N _i( E ), respectively.

In the virtual-voltage method, both T and $\{{\phi}_{i}\}_{i=∼1}^{m}$ are determined iteratively by means of the following two steps:

The linear boundary-value problem (7)–(10) is solved for δT and $\{{\delta}{\phi}_{i}\}_{i=∼1}^{m}$ .

Both T and 𝜙_i (i = 1,2, …, m) are updated by $\begin{eqnarray}T\leftarrow T+{\delta}T,\quad {\phi}_{i}\leftarrow {\phi}_{i}+{\delta}{\phi}_{i}\quad (i=1,2,\ldots ,m).\end{eqnarray}$

The above two steps are repeated until both convergence of T and that of $\{{\phi}_{i}\}_{i=∼1}^{m}$ are achieved.

4. Linear-system solver

4.1. Linear system at each iteration cycle of Newton method

If the linear boundary-value problem (7)–(10) is discretized with the finite element method with n nodes, the following linear system is obtained: $\begin{eqnarray}\displaystyle \left[\begin{array}{@{}cc@{}}A(\boldsymbol{T}) & C\\ D^{T}(\boldsymbol{T}) & O\end{array}\right]\left[\begin{array}{@{}c@{}}{\delta}\boldsymbol{T}\\ \boldsymbol{\it\lambda}\end{array}\right]=\left[\begin{array}{@{}c@{}}\boldsymbol{G}(\boldsymbol{T},\boldsymbol{{\phi}})\\ \boldsymbol{ c }(\boldsymbol{T })\end{array}\right]. & & \displaystyle\end{eqnarray}$ (11) Here, T and δ T are nodal vectors corresponding to T and δT, respectively, and 𝜙 is defined as $\boldsymbol{{\phi}}=\left(\begin{array}{@{}cccc@{}}{\phi}_{1} & {\phi}_{2} & \cdots ∼ & {\phi}_{m}\end{array}\right)^{T}$ . In addition, 𝝀 is a (k + m)-dimensional unknown vector whose ith component is δ𝜙_i (i = k +1, k +2, …, k + m). Here, 2k denotes the number of elements adjacent to crack surfaces and it is assumed to satisfy n ≫ k. In the following, k is called the number of crack elements. Also, $A(\boldsymbol{T})\in \mathbb{R}^{n\times n}$ , $C\in \mathbb{R}^{n\times (k+m)}$ are derived from the weak form equivalent to (7)–(9), whereas D ^T( T ) δ T = c ( T ) is a discretized form of (10) and the essential boundary condtion, $\partial$ T∕ $\partial$ s = 0 on C ₁, C ₂, …, C _m. Furthermore, $\boldsymbol{G}(\boldsymbol{T},\boldsymbol{{\phi}})\in \mathbb{R}^{n}$ and $\boldsymbol{c}(\boldsymbol{T})\in \mathbb{R}^{(k+m)}$ are vector-valued functions.

The linear system (11) changes depending both on an iteration cycle of the Newton method and on time. Hence, throughout this section, two linear-system solvers are applied to (11) for the 1st iteration cycle at t = Aw∕(1200v) and their performances are compared with each other. Numerical computations were carried out on FUJITSU PRIMEHPC FX100 of LHD Numerical Analysis Server in National Institute of Fusion Science.

4.2. Conventional method

Since both n ≫ k and n ≫ m are satisfied, most elements of the coefficient matrix in (11) are occupied by the dense submatrix A ( T ). Hence, if the LU decomposition is applied to the solution of (11), its computational cost amounts to O (n ³). In other words, the LU decomposition is extremely time-consuming and, hence, it is inappropriate for the numerical solution of (11). For this reason, the authors have applied GMRES to the solution of (11). In the present paper, the numerical solution of (11) by GMRES is called the conventional method.

Let us first investigate residual histories of the conventional method. To this end, residual histories are determined for various values of k and are depicted in Fig. 2. This figure indicates that the convergence property of GMRES is degraded with an increase in k. For the purpose of quantitatively investigating this tendency, we use the number of iterations required for convergence of GMRES. In the following, it is called the convergent iteration number, N _c.

Fig. 2.

Residual histories of the conventional method for the case with n = 11649.

Next, the dependence of the convergent iteration number on k is numerically determined and is shown in Fig. 3. We see from this figure that, for the conventional method, the convergent iteration number increases linearly with k. Hence, even for the case where the same node distribution is used, an increase in the number k of crack elements will raise the CPU time required for the solution of (11).

Fig. 3.

Dependences of the convergent iteration number N _c on the number k of crack elements for the case with n = 11649.

4.3. Variable-reduction method

As mentioned above, the computational cost of the conventional method strongly depends on the number k of crack elements. Such a strong dependence of the computational cost on k is attributable to the fact that 𝝀 is contained in (11).

Let us decompose (11) into two equations so that δ T and 𝝀 can be separately determined. For this purpose, the QR factorizations [10] of C and D ( T ) are first computed as follows: $\begin{eqnarray}C=Q_{C}R_{C}P_{C}^{T},\quad D(\boldsymbol{T})=Q_{D}(\boldsymbol{T})R_{D}(\boldsymbol{T})P_{D}^{T}(\boldsymbol{T}),\end{eqnarray}$ where $Q_{C},Q_{D}(\boldsymbol{T})\in \mathbb{R}^{n\times (k+m)}$ are matrices such that $Q_{C}^{T}Q_{C}=I$ and Q _D ^T( T )Q _D( T ) = I. Here, I denotes an identity matrix. Furthermore, $R_{C},R_{D}(\boldsymbol{T})\in \mathbb{R}^{(k+m)\times (k+m)}$ and $P_{C},P_{D}(\boldsymbol{T})\in \mathbb{R}^{(k+m)\times (k+m)}$ denote upper triangular matrices and permutation matrices, respectively.

By using projection matrices, $U_{C}\equiv I-Q_{C}Q_{C}^{T}$ , F ( T ) ≡ Q _C Q _D ^T( T ) and U ( T ) ≡ I − F ( T ), we can easily prove that the linear system (11) is equivalent to the following two equations: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle A^{\ast }(\boldsymbol{T})∼{\delta}\boldsymbol{T}=\boldsymbol{G}^{\ast }(\boldsymbol{T},\boldsymbol{{\phi}}),\end{eqnarray}$ (12) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle \boldsymbol{\it\lambda}=P_{C}R_{C}^{-1}Q_{C}^{T}[\boldsymbol{G}(\boldsymbol{T},\boldsymbol{{\phi}})-A(\boldsymbol{T})∼{\delta}\boldsymbol{T}].\end{eqnarray}$ (13) Here, the matrix $A^{\ast }(\boldsymbol{T})\in \mathbb{R}^{n\times n}$ and the vector $\boldsymbol{G}^{\ast }(\boldsymbol{T},\boldsymbol{{\phi}})\in \mathbb{R}^{n}$ are given by $\begin{eqnarray}A^{\ast }(\boldsymbol{T})=U_{C}^{T}∼A(\boldsymbol{T})∼U(\boldsymbol{T})+F(\boldsymbol{T}),\quad \boldsymbol{G}^{\ast }(\boldsymbol{T},\boldsymbol{{\phi}})=U_{C}^{T}[\boldsymbol{G}(\boldsymbol{T},\boldsymbol{{\phi}})-A(\boldsymbol{T})∼\boldsymbol{c}^{\ast }(\boldsymbol{T})]+\boldsymbol{c}^{\ast }(\boldsymbol{T}),\end{eqnarray}$ where c ^∗( T ) = Q _C R _D ^−T( T )P _D ^T( T ) c ( T ). Note that 𝝀 is not contained in (12).

By solving (12) for δ T and, subsequently, substituting δ T into (13), 𝝀 is also calculated. In the actual calculation, (12) is numerically solved with GMRES. Throughout the present study, this method is called the variable-reduction method.

We first investigate the convergence property of the variable-reduction method. To this end, the convergent iteration number is numerically determined as a function of k and is also plotted in Fig. 3. For the variable-reduction method, the convergent iteration number is little affected by k.

Fig. 4.

Dependences of the CPU time τ on the number n of nodes. Here, (a) L _c = 15 mm and (a) L _c = 45 mm.

Next, we investigate the speed of the variable-reduction method. The CPU times required for solving (11) either with the conventional method or with the variable-reduction method are measured, and they are depicted in Figs 4(a) and 4(b). These figures indicate that the variable-reduction method is always faster than the conventional method. In other words, the variable-reduction method shows the speedup effect as compared with the conventional method. Moreover, the speedup effect becomes remarkable with an increase either in the number of nodes or in a crack size.

5. Conclusion

We have developed a novel method for solving a linear system in the shielding current analysis of a cracked HTS film. After discretized with respect to time and space, the initial-boundary-value problem of the shielding current density reduces to a linear system at each iteration cycle of the Newton method. In order to solve the linear system efficiently, the variable-reduction method is proposed. Furthermore, its performance is numerically investigated. Conclusions obtained in the present study are summarized as follows:

–

Although the convergence property of the conventional method is significantly degraded with increasing crack size, that of the variable-reduction method is hardly influenced by a crack size.

–

The speed of the variable-reduction method is faster than that of the conventional method. This speedup effect becomes remarkable with an increase either in the number of nodes or in a crack size.

Footnotes

Acknowledgements

This work was supported in part by Japan Society for the Promotion of Science under a Grant-in-Aid for Scientific Research (C) No. 15K05926. A part of this work was also carried out with the support and under the auspices of the NIFS Collaboration Research program (NIFS17KNTS051, NIFS16KNXN341).

References

Yoshida

Uesaka

and Miya

, Magnetic-field and force analysis of high-T _c superconductor with flux-flow and creep, IEEE. Trans. Magn. 30 (1994), 3503–3506.

Kamitani

Takayama

and Ikuno

, High-speed method for analyzing shielding current density in high-temperature superconductor, IEEE. Trans. Magn. 47 (2011), 1138–1141.

Kameni

Boubekeur

Alloui

Bouillault

Lambretchs

and Geuzaine

, A 3-D semi-implicit method for computing the current density in bulk superconductors, IEEE. Trans. Magn. 50 (2014), 7009204(4pp).

Kamitani

Takayama

and Ikuno

, Virtual voltage method for analyzing shielding current density in high-temperature superconducting film with cracks/holes, IEEE Trans. Magn. 49 (2013), 1877–1880.

Kamitani

Takayama

and Saitoh

, Numerical simulation of permanent magnet method: Applicability to crack detection in high-temperature superconducting film, IEEE Trans. Magn. 51 (2015), 7000504(4pp).

Kamitani

Takayama

and Saitoh

, High-speed shielding current analysis in high-temperature superconducting film with cracks, IEEE. Trans. Magn. 52 (2016), 7202404(4pp).

Makong

Kameni

Masson

Lambrechts

and Bouillault

, 3-D modeling of heterogeneous and anisotropic superconducting media, IEEE. Trans. Magn. 52 (2016), 7205404(4pp).

Dianati

Heydari

and Afsari

S.A.

, Analytical computation of air-gap magnetic field in a viable superconductive magnetic gear, IEEE Trans. Appl. Supercond. 26 (2016), 5205612(4pp).

Hattori

Saito

Takano

Suzuki

Yamada

Takayama

Kamitani

and Ohshima

, Detection of smaller J _c region and damage in YBCO coated conductors by using permanent magnet method, Physica C 471 (2011), 1033–1035.

10.

Golub

G.H.

and Van Loan

C.F.

, Matrix Computations, 4th edn, Johns Hopkins University Press, MD, USA, 2013, pp. 246–259.