Application of parallel processing technique to shielding current analysis on HTS thin film

Physica C 463–465 (2007) 1013–1016 www.elsevier.com/locate/physc

Application of parallel processing technique to shielding current analysis on HTS thin film Soichiro Ikuno a

a,*

, Teruou Takayama b, Atsushi Kamitani

b

School of Computer Science, Tokyo University of Technology, Katakura 1404-1, Hachioji, Tokyo 192-0982, Japan b Faculty of Engineering, Yamagata University, Johnan 4-3-16, Yonezawa, Yamagata 992-8510, Japan Received 1 November 2006; received in revised form 1 December 2006; accepted 19 March 2007 Available online 26 May 2007

Abstract Parallelized element-free Galerkin (EFG) method using OpenMP is evaluated and the evaluation of AC loss of axisymmetric high-Tc superconductors (HTSs) by using EFG is presented. Formulation of the electromagnetic behavior of shielding current density in an HTS gives a system of time-dependent integro-differential equations. After discretizing the weak form of the problem based on EFG, the axisymmetric simulation code for analyzing the time evolution of the shielding current density on HTS is developed. In the present study, the simulation code for analyzing the time evolution of the shielding current density in HTS by using EFG is developed. Results of the computation shows that the CPU time of parallelized EFG in the case of 4PU is 3.8 times faster than that of 1PU. Furthermore, the AC loss is calculated numerically by use of the values of the shielding current density.  2007 Elsevier B.V. All rights reserved. PACS: 07.05.Tp; 41.20.Gz; 74.70.b; 85.25.j Keywords: Magnetic shielding; AC loss; YBCO; Numerical simulation; Galerkin method

1. Introduction The time dependent spatial distribution of the shielding current density in high-Tc superconductors (HTSs) is necessary information to calculate the AC loss, levitation force and other physical quantities in the HTS numerically [1]. Various methods for calculating shielding current density have been proposed [2,3]. Formulation of the electromagnetic behavior of the shielding current density in an HTS gives a system of time-dependent integro-differential equations. The system is solved with the J–E constitutive relation, which is dependent on applied magnetic flux. The behavior can be determined by solving the initial-boundary-value problem of the system. A mesh generation procedure is necessary for using the general finite element method (FEM). Generally, it costs *

Corresponding author. Tel./fax: +81 426 37 2793. E-mail address: [email protected] (S. Ikuno).

0921-4534/$ - see front matter  2007 Elsevier B.V. All rights reserved. doi:10.1016/j.physc.2007.03.471

huge amount of time to divide the region into a set of finite elements. On the other hand, the element-free Galerkin (EFG) method [4] does not require finite elements of a geometrical structure. However, it takes tremendous CPU time to generate the coefficient matrix of linear system instead of unnecessary mesh generation procedure. The purpose of the present study is to develop numerical code for analyzing shielding current density in HTSs. In the present study, the power law is employed to construct the J–E constitutive relation; the independent variable in the power law depends on applied magnetic flux. Moreover, the parallelized EFG is evaluated and AC losses against the applied magnetic flux whose frequencies are from 100 Hz to 500 Hz are calculated. 2. Governing equation and discretization Let us first explain the governing equation of shielding current density in a melt-powder-melt-growth processed

1014

S. Ikuno et al. / Physica C 463–465 (2007) 1013–1016

YBCO (MPMG-YBCO) superconductor. The HTS is placed in homogeneous AC magnetic field with its flux parallel to the thickness direction. For simplicity, we assume that the HTS is disk-shaped whose radius is R and whose thickness is D, and the area of the circular cross-section is constant along the thickness direction. Thus, we can treat the problem as axisymmetric. It is known that an MPMG-YBCO superconductor has crystallographic anisotropy in its critical current density jc [5]. By taking this fact into account, we apply multiplethin-layer approximation to the problem; we also assume the shielding current density does not flow across any interface of two adjacent layers. Under the conditions above, the HTS disk is regarded as being composed of M pieces of thin layers of thickness 2e. In other words, the thickness of the HTS, D, satisfies the equation D = 2e · M. Throughout this paper, let us use the cylindrical coordinate (r, h, z) by taking the symmetry axis of the HTS as zdirection. In terms of the coordinate, the applied magnetic flux density B0 is written as B0 = B0sin 2pftez where B0 denotes amplitude of flux density and f denotes frequency of flux and ez denotes the unit vector in the z-direction. Under these assumptions, there exists a scalar function Sp(r, t) for the pth layer such that the shielding current density in this layer, jp, satisfies the equation j p ¼ 1e rS p  ez . Here, the subscript p indicates values in the pth layer. Using this scalar function, the behavior of the shielding current density can be expressed by the following timedependent integro-differential equation: M Z R X oS q 0 l0 oS p l0 dr þ r0 Qpq ðr; r0 Þ ot e ot 0 q¼1 ¼

oðB 0  ez Þ  ðr  E p Þ  ez ; ot

ð1Þ

where l0 denotes the magnetic permeability of vacuum. The function Qpq(r, r 0 ) is defined as Qpq ðr; r0 Þ ¼ 

1 X 1 X 1 mþn mn pffiffiffiffiffi ð1Þ k mn pq Kðk pq Þ; 2 0 4pe rr m¼0 n¼0

ð2Þ

where K(x) denotes x’s complete elliptic integral of the first kind, and k mn pq is defined by 2 ðk mn pq Þ ¼

4rr0 ðr þ

r0 Þ2 fðz

p

 zq Þ þ ½ð1Þm  ð1Þn eg2

;

ð3Þ

where zp denotes the z coordinate of the middle of the pth layer. For the initial and boundary conditions to Eq. (1), we assume Sp(r, 0) = Sp(R, t) = 0. The initial-boundary-value problem of Eq. (1) is solved with the J–E constitutive relation which is defined as E p ¼ Eðjj p j; BÞ

jp ; jj p j

ð4Þ

where B denotes the magnitude of the applied magnetic flux density. In the present study, the function E(j, B) is given by the power law and it can be written in the form

 Eðj; BÞ ¼ Ec

j jc ðBÞ

N ;

ð5Þ

where Ec denotes the critical electric field. In this paper, N is fixed as N = 19. The B-dependence of jc is also defined as jc ðBÞ ¼ jc0 ½BK =ðB þ BK Þ;

ð6Þ

where jc0 and BK are constants. The critical current density is a function of B/BK and the maximum magnitude of applied magnetic flux density is B0. That is to say, the Bdependence of critical current density jc is characterized by the parameter B0/BK. By applying the backward Euler method to Eq. (1), the system is discretized with respect to time, and it is transformed to the boundary-value problem. The problem is expressed as a weak form, and it is discretized with respect to space by using the EFG method. The obtained nonlinear system is written as G(s, k) = 0. Here the nodal vectors s and k correspond to the scalar function Sp, the electric field and the Lagrange multiplier, respectively. The nodal vector s can be determined by solving the nonlinear system G(s, k) = 0. As we mentioned above, it takes tremendous CPU time to generate the coefficient matrix of the system using EFG. That is to say, the speedup of this procedure leads to the reduction of the experiments time. Note that B-dependent critical density and the power law are adopted for the J–E constitutive relation. Thus, we use the adaptively deaccelerated Newton method (ADNM) [6] for the nonlinear system. 3. Evaluation of parallelized EFG Firstly, let us investigate the ratio of CPU time that the procedure in EFG spends. For the simplicity, Poisson problem is adopted for the evaluation. We investigate the ratio of CPU time that each procedure spends by using profiling mode of the compiler. Note that the LU decomposition method is employed for the solver for the linear system. As a result, it has been understood that CPU time used to procedure for making the coefficient matrix is about 40%. In addition, 50% is used for calculating the linear system by use of LU decomposition method, and 10% is used for I/O and other procedure. From the above result, the speedup of procedure for making the coefficient matrix leads to the reduction of CPU time. Thus, let us parallelize the procedure for making the coefficient matrix using the OpenMP. The OpenMP is a one of application programming interface (API) for multiprocessing programming in C/C++ and Fortran on shared memory multiprocessor architecture. The OpenMP consists of a set of compiler directives, library routines and environment variables. Fig. 1 shows the CPU time for making the coefficient matrix as function of total number of node N that is used for EFG. The evaluation environment are shown as follows; machine: AMD Opteron 270 (Dual Core, 2 GHz) · 2 = 4 processor, memory: 4GB, OS: FreeBSD 6.1-STABLE,

S. Ikuno et al. / Physica C 463–465 (2007) 1013–1016

Fig. 1. The CPU time for making the coefficient matrix as function of total number of node N in the case of Nc = 1024, Ng = 32. Here Nc denotes number of cell that divides the region and Ng denotes number of integral point for Gauss–Legendre numerical integration.

1015

D = 2 mm, M = 6, T = 77 K, jc = 1.3 · 107 A/m2 and Ec = 0.1 mV/m. First, we show the time evolution of the shielding current density. The shielding current density of h component is plotted as functions of 2ft in Fig. 3. The applied magnetic flux is also plotted in the same figure. Here the value of shielding current density is calculated at r/R = 0.6 in the first layer. We can see some lagging phases between the applied magnetic flux and the shielding current density, but the value of shielding current density changes periodically and recurrently as the applied magnetic flux density changes. Next, we investigate the AC losses of the HTS. The total time-dependent AC loss is calculated by the following equation: M Z R X Qtime ¼ 4pe rðE p  j p Þdr: ð7Þ p¼1

0

Note that values of the physical quantities are evaluated in the first layer of the HTS through this paper. In Fig. 4, we show the time dependence of total AC loss. The value of AC loss changes periodically and recurrently in doublephase as the shielding current density changes. Note that

Fig. 2. The speedup ratio as function of number of processors in the case of N = 4096, Nc = 1024, Ng = 32. Here Nc denotes number of cell that divides the region and Ng denotes number of integral point for Gauss– Legendre numerical integration.

compiler: GNU fortran compiler 4.2.0, compiler option: AO2–fopenmp. Fig. 1 shows that the CPU time for making the coefficient matrix increases as the total number of nodes increases in both cases. Moreover, the CPU time rapidly increases at large N in the case of 1PU. On the other hand, the exponential growth of the CPU time is suppressed in the case of 4PU. The speedup ratio is plotted as function of number of processors in Fig. 2. We see form this figure that, for all the cases, good performances are observed. 4PU demonstrates 3.8 times faster than 1PU. From this result, it can be concluded that there is an enough parallelized effect.

Fig. 3. Time evolution of the z component of the applied magnetic flux B0z (dotted line) and the shielding current density (solid line) in case of f = 300 Hz, B0 = 0.3 T.

4. AC loss analysis In this section, we apply the method given above to solve the initial-boundary-value problem to determine the behavior of the shielding current density in the HTS. Moreover, we calculate the AC loss from the shielding current density. Throughout this paper, the geometrical and the physical parameters are fixed as follows: R = 20 mm,

Fig. 4. Total time dependent AC loss against nondimensional time 2ft in case of f = 300 Hz, B0 = 0.3 T.

1016

S. Ikuno et al. / Physica C 463–465 (2007) 1013–1016

increases as f increases as well as Fig. 6. However, compared to the values in Fig. 5, Q increases in a linear fashion in Fig. 6. This result indicates that the frequency of magnetic flux density affects directly to AC loss compared to its magnitude. 5. Conclusions

Fig. 5. Time-averaged total AC loss against the value of B0. (A) f = 100 Hz, (B) f = 300 Hz, (C) f = 500 Hz.

Fig. 6. Time-averaged total AC loss against the frequency of the applied magnetic flux density. (A) B0 = 0.1 T, (B) B0 = 0.3 T, (C) B0 = 0.5 T.

there is not much difference between the values of AC loss in any two layers. Finally, we investigate the influence of the magnitude and frequency of the applied magnetic flux density. In order to discuss it, time-averaged total AC loss, Q, is calculated using the following equation Z 1 T Q¼ Qtime dt: ð8Þ T 0 Q is plotted as function of B0 and f in Figs. 5 and 6. Fig. 6 shows that Q increases smoothly as B0 increases. This figure also indicates that the increase of f leads to increase of Q. In order to investigate the influence of f to Q in more depth, Q is plotted as function of f in Fig. 6. In this figure Q

We have developed the numerical code for analyzing the time evolution of the shielding current density in axisymmetric HTSs using the EFG method and evaluated the AC loss against the value of B0 and frequency of applied magnetic flux density. Conclusions obtained in the present study are summarized as follows: (1) The CPU time rapidly increases at large N in the case of 1PU. On the other hand, the exponential growth of the CPU time is suppressed in the case of 4PU. The speedup ratio of 4PU is 3.8 times larger than that of 1PU. (2) The lagging phases occur between the applied magnetic flux and the shielding current density, but the values of shielding current changes periodically and recurrently as the applied magnetic flux density changes. (3) The value of AC loss changes periodically and recurrently in double-phase as the shielding current density changes. (4) Increase of the value of the frequency and the value of the magnetic flux density lead to increase of the AC loss. (5) The frequency of magnetic flux density affects directly to AC loss compared to its magnitude.

References [1] N. Enomoto, T. Izumi, N. Amemiya, IEEE Trans. Appl. Supercond. 15 (2005) 1574. [2] Y. Yoshida, M. Uesaka, K. Miya, Int. J. Appl. Electromagn. Mater. 5 (1994) 83. [3] M. Uesaka, A. Suzuki, N. Tanaka, Y. Yoshida, K. Miya, Cryogenics 35 (1995) 243. [4] T. Belytschko, Y.Y. Lu, L. Gu, Int. J. Numer. Meth. Eng. 37 (1994) 229. [5] M. Murakami, S. Gotoh, H. Fujimoto, K. Yamaguchi, N. Koshizuka, S. Tanaka, Supercond. Sci. Technol. 4 (1991) S43. [6] A. Kamitani, S. Ikuno, T. Yokono, IEICE Trans. Electron. E 87-C (2004) 101.