Fast solution of the superconducting dynamo benchmark problem
11 Fast solution of the superconducting dynamo benchmark problem
Leonid Prigozhin and Vladimir Sokolovsky J. Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Sde Boqer Campus 84990, Israel Physics Department, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel E-mail: [email protected] and [email protected],
Abstract
A model of high temperature superconducting dynamo, a promising type of flux pumps capable of wireless injection of a large DC current into a superconducting circuit, has recently been chosen as an applied superconductivity benchmark problem and solved using eight different finite element methods (Ainslie et al
Supercond. Sci. Technol. Keywords: superconducting dynamo, coated conductor, numerical solution, Chebyshev polynomials, singular integral equation Introduction
The high- c T superconducting (HTS) dynamos are devices capable of inducing a large DC current in a superconducting circuit without the cryogenic losses associated with non-superconducting current leads [1-5]. Mathematical models, explaining the DC voltage in an HTS strip (the dynamo stator) generated by the dynamo-rotor-attached permanent magnet (PM), have been proposed recently, see [6-10] and the references therein. Available mathematical descriptions are either simplified circuit-type models of the HTS dynamo (e.g. [2, 4, 6]) or consider an infinite superconducting strip subjected to the time-varying non-uniform field of a rotating long PM [7-10]. Models of the latter type are less simplified and better reproduce the non-trivial features of the HTS dynamos. Such HTS dynamo model, with a coated conductor stator characterized by the power current-voltage law with a field-independent critical current density, has been proposed as a new benchmark problem in applied superconductivity and solved numerically by eight different finite element methods in [11]. These authors computed the generated open-circuit voltage, compared the computation times needed to the methods employed, and made their simulation results available as a supplementary material. Here we show that an accurate numerical solution of this benchmark problem can be obtained using the method of lines for integration in time and Chebyshev polynomials for the approximation in space. For HTS strips, strip stacks, and pancake coils a similar approach has been proposed in [12]. As in [12], our choice of Chebyshev polynomials for spacial approximation is based on two factors: remarkably fast convergence of the interpolation by Chebyshev polynomial expansions [13] and the ease of using such expansions for solving integral equations with a singular Cauchy-type kernel (see, e.g., [14, 15] ). The derived method is much faster than all methods in [11]. Numerical simulations in our work were done using Matlab R2020a on a PC with the Intel (R) Core (TM) i7-9700 CPU@ 3,00 Hz and 32 GB RAM. The benchmark problem
Following [11], we assume that an infinitely long PM with the remanent flux density r B rotates counter clockwise past a thin stationary HTS strip in the open-circuit configuration (Figure 1); all HTS dynamo parameters are listed in Table 1. Figure 1. A scheme of an HTS dynamo: the geometry of the benchmark problem.
Table 1. HTS dynamo benchmark parameters (as in [11]). permanent magnet (PM) width, w v
12 mm active length, l r B a
12 mm thickness 1 μm critical current, c I
283 A power value, n
20 rotor external radius, R
35 mm air gap, d f If a long permanent magnet is in the initial position (Figure 1), its magnetic field ( , )
Ty z
H H H at a point ( , ) T y z p outside the magnet is, see [16], ( ) ( ) ( ) ( )( ) ln ln ,4 ( ) ( ) ( ) ( )2 ( ) 2 ( )( ) arctan arctan ,2 ( ) ( ) yz B y w z v y w z vH y w z v y w z vB v y w v y wH y w z v y w z v pp where is the permeability of vacuum, z z R v , R is the rotor radius, v and w are the PM height and width, respectively. After rotation of the rotor by an angle the field at the point p becomes PM 0 ( , ) ( ) ( ( ) ),
A A
H p H p (1) where cos sin( ) sin cos A is the matrix of rotation. Neglecting the strip thickness, we present the strip as s , , x a y a z z , where a is the strip width, s z R d , and d is the air gap. It is assumed the parallel-to-strip electric field component E and the strip sheet current density J are directed along the strip (parallel to the x -axis) and obey the power law,
10 c c . n J JE E J J (2) Here, as in [11], the field-independent sheet critical current density is c c / 2
J I a , where c I is the strip critical current, E V/m, the power n is a constant. By the Faraday law, at the strip points s s ( , ) T y z p we have z HEy t . (3) Here the normal-to-strip component z H of the total magnetic field is a superposition of the PM and the strip current field components: PM, s az z a
J t uH t y H t uy u p . (4) Since for all strip points s z z , to simplify our notations we will omit the dependence on z and write PM, ( ( ), ) z H t y instead of PM, s ( ( ), ) z H t p . Finally, if a source of current is present, the additional condition is ( , )d ( ) aa J t y y I t , (5) where ( ) I t is a given function. The strip open-circuit condition in the benchmark problem means the transport current I is zero. Jointly with the initial condition (0, ) 0 J y , equations (1)-(5) fully describe the HTS dynamo benchmark problem [11]. Experimentally, the instantaneous open-circuit voltage at the ends of a stator strip segment has been measured in [2, 7, 9, 10]. Taken for a segment containing the close-to-rotor active zone of the stator, this voltage and its time averaged DC value are the most important characteristics of an HTS dynamo. In the benchmark model the voltage is suggested to calculate as ( ) ( , )d ,2 aa lV t E t y ya (6) where l is the length of the active zone (Table 1). Numerical scheme
Rescaling the variables, c c 000 c c / , ( , ) ( , ) / , / ,( , ) ( , ) / , , ,2 y z y z
J J J H H H H J E E EE Iy z y z a t t IJ a aJ omitting the primes, and differentiating equations (4) and (5) in respect to time, we obtain in dimensionless form for
0, 1 1 t y : | | , n E J J (7) z J t u Eu Hy u y (8) ( , )d 2 ( ). J t y y I t (9) Here the dot above variables means time derivative, the singular integral is understood in the principal value sense. As in [12], we seek an approximation to the sheet current density in the form of a weighted expansion in Chebyshev polynomials of the first kind, ( ) ( )( , ) 1
N i ii t T yJ t y y , (10) and derive a system of ordinary differential equations (ODE) for N unknowns, the values ( ) ( , ) k k J t J t y in N mesh points, the roots of the Chebyshev polynomial ( ) cos arccos( ) N T y N y : k y k k NN . (11) Using such meshes suppresses the Runge phenomenon that makes polynomial interpolation on the uniform meshes unstable. Substituting (10) into equation (8) we use the well known relation [17] d1 m m T u u U yy u u , (12) where [ 1,1] y , ( ) sin( arccos ) / sin(arccos ) m U y m y y for
1, 2,... m are Chebyshev polynomials of the second kind, and U . This yields, for
1, 2,..., k N , k N i i k zi y
EU y Hy . (13) Substituting (10) also into equation (9) and noting that m mT y y my , we find I t . Let the k J values be known at time . t To find the right hand side of equations (13) we calculate the electric field (7) at the mesh nodes, | | nk k k E J J , and, approximately, find / E y at the same points by differentiating the Chebyshev expansion interpolating the k E values (see Appendix A). Furthermore, PM, PM,PM, d 2d z zz
H HH ft . Computing
PM, / z H for ( ) t at the mesh points k y (Appendix B), we calculate PM, ( ) / k k z y y t E y H . (14) It is not difficult to construct the expansion in Chebyshev polynomials of the second kind, ( ) N m mm
U y , such that ( ) k k y (Appendix A). Equations (13) yield then i i i N . Finally, we set ( ) ( )( ) 1
N i i kik k t T yJ t y . (15) The ODE system (15) should be integrated in time. In our work we employed a standard Matlab ODE solver, ode15s , with the default parameters. Solution of the benchmark problem
To compare our numerical simulation results to those in [11] we return to dimensional variables. As in [11], we assumed zero transport current, simulated ten full rotor rotations, and presented the ( )
V t curves for the second cycle to avoid the transient dynamics during the first cycle.
First, to estimate the accuracy and convergence rate of our scheme, we solved the problem on different meshes and found that, visually, solutions for N and N are very close (Figure 2). Comparison with the latter, our most accurate solution was used to estimate the relative error in L norm for several meshes (Table 2). Figure 2. Dynamo generated voltage. Numerical solutions for 100 and 1000 mesh nodes. Table 2. Simulation results. N Time per cycle (s) Estimated error (%) 100 1.2 3 200 4.8 1 400 33 0.3 1000 612 -- No accuracy estimates for other methods are available and, although the presented ( )
V t curves ([11], Figure 3) look similar, small deviations can be noticed in some of them. Using the provided supplementary data we compared our solution to those of the four fastest finite element methods in [11] called there MEMEP, SEG-H, VIE, and H-A (Figure 3). The reported computation time for the variational MEMEP method is less than for other methods (about 15 seconds per cycle). However, the presented MEMEP solution is not very accurate and differs from ours and other solutions in [11]. Significantly more time per cycle was needed in [8], where the same method was used with finer grids and, probably, more accurate solutions obtained. Each of the three other methods needed more than 1 minute per cycle and their solutions are close to ours. Our method is, however, much faster.
Figure 3.
The benchmark problem. Solution by the Chebyshev polynomial-based method versus solutions by the four fastest finite element methods from [11]. Extension of the benchmark model
Our method is easily extended to the model with a field-dependent sheet critical current density, e.g., c0c 1 2 20 ,1 y z JJ h H H (16) where c0 0 , , J h are fitting parameters. The parallel-to-film magnetic field component y H , discontinuous on an infinitely thin film, is usually replaced by the corresponding component of the external field ( PM, y H in our case). To compute z H on each time step, the mesh values k k J y are interpolated by the Chebyshev expansion ( ) N k kk
T y (Appendix A). This presents J in the form (10); substituting (10) into equation (4) and using (12) one obtains
1( , ) ( ( ), ) ( )2
Nz k z k i i ki
H t y H t y U y . Finally, the mesh values of c J are calculated using (16) and the electric field values k E are determined by the dimensionless form of equation (2). In our simulations (Figure 4) we chose c0 J equal to the constant c J from the benchmark problem, , and three different values for h . For h we have c c0 J J and, the smaller h is, the stronger is the dependence of the sheet critical current density on the magnetic field and the higher is the time-averaged DC voltage. Qualitatively, the voltage curves in Figure 4 are similar to the experimental ones, see [2, 7, 10]. A closed-circuit problem with an HTS dynamo connected to an electronically controlled current supply has been studied, both experimentally and numerically, in [10]. Numerical simulation in this work employed the H -formulation, and so a two-dimensional finite element mesh was needed in the plane orthogonal to the strip; the strip thickness has had to be artificially increased from 1 to 100 μm . A one-dimensional formulation, like the one used in our work, is expected to serve a basis for a more efficient numerical scheme for such problems. We computed the instantaneous HTS dynamo voltage for several transport currents assuming the field dependent critical sheet current density (16) with c h J . In these examples the transport current changed linearly with time during the first quarter of the first cycle and then remained constant; numerical solutions are shown in Figure 5. For the examples in this Section the computation for N took about 6.5 second per cycle and, as above, the solutions are shown for the second cycle. Figure 4. Simulation results for models with field-dependent critical sheet current density. Figure 5. Simulation results for problems with transport current. Here the c ( ) J H dependence (16) is employed with c h J . Discussion
Devices for wireless injection of DC supercurrents into closed superconducting circuits open a way to HTS magnets without the thermally inefficient metal current leads and attracted much recent interest (see the review [18]). The HTS dynamo considered in our work is a device of this type. Its physical mechanism, leading to creation of an open-circuit DC output voltage has been explained in [2, 7-10] and a corresponding mathematical model was proposed [7-10]. The model is able to describe, at least qualitatively, the main features of this device. Numerical simulations based on this model were used in these works to investigate the influence of various factors (width and critical current of the coated conductor stator, air gap, PM rotation frequency, etc.) on the HTS dynamo performance. In [11], this model problem was chosen as a benchmark problem and solved by several popular finite element methods; a comparison of the efficiency of these methods has been presented. Although useful, such a comparison does not fully characterize these methods, especially those derived for more general (multidimensional) problems. However, as shown in our work, for this specific benchmark and similar problems the Chebyshev polynomial-based method is more than an order of magnitude faster than all finite element methods in [11]. Partly, this is related to the efficient treatment of the singular Cauchy integral term in equation (4), based on the relation (12). Here we solved the HTS dynamo open circuit benchmark problem as well as the closed - circuit problems with a current source and problems with field - dependent sheet critical current density. In the latter case we used the most often employed form (16) of this dependence. However, the dependence can be easily replaced by, e.g., interpolation using the experimentally obtained data, as in [7, 8, 10]. The method can be used also to model the traveling wave flux pumps [18] and multi-PM dynamos [5, 18]. ORCID iD
Leonid Prigozhin https://orcid.org/0000-0002-2448-4471
Vladimir Sokolovsky https://orcid.org/0000-0003-4887-413X
References [1] Hoffmann C, Pooke D, and Caplin A D 2010 Flux pump for HTS magnets
IEEE Trans. Appl. Supercond. T c superconducting dynamo Appl. Phys. Lett.
108 122601 [3] Hamilton K, Mataira R, Geng J, Bumby C, Carnegie D and Badcock R 2020 Practical estimation of HTS dynamo losses
IEEE Trans. Appl. Supercond. et al. Rotating permanent magnets based flux pump for HTS no-insulation coil
IEEE Trans. Appl. Supercond. IEEE Trans. Appl. Supercond. Supercond. Sci. Technol. T c superconducting dynamo Appl. Phys. Lett.
Supercond. Sci. Technol. T c superconducting dynamos IEEE Trans. Appl. Supercond. T c superconducting dynamo: Models and experiment Phys. Rev. Appl . , et al. T c superconducting dynamo Supercond. Sci. Technol. Supercond. Sci. Technol. Approximation theory and approximation practice (Philadelphia, PA: Siam) [14] Eshkuvatov Z, Long N N and Abdulkawi M 2009 Approximate solution of singular integral equations of the first kind with Cauchy kernel
Appl. Math. Lett . Comm. Adv. Math. Sci. Permanent magnet and electromechanical devices: materials, analysis, and applications (San Diego, CA: Academic Press) [17] Mason J C and Handscomb D C 2002
Chebyshev polynomials (Baca Raton, FL: CRC press) [18] Zhai Y, Tan Z, Liu X, Shen B, Coombs T A and Wang F 2020 Research Progress of Contactless Magnetization Technology: HTS Flux Pumps
IEEE Trans. Appl. Supercond. Appendix A. Operations with Chebyshev interpolating expansions
Our numerical scheme makes use of the following linear transformations: ( i ) from function values at the Chebyshev mesh (11) to coefficients of the interpolating Chebyshev polynomial expansion and vice versa; ( ii ) from the mesh values of the electric field E to the approximate values of / E y at the same mesh computed as derivatives of the interpolating polynomial; ( iii ) from the mesh values of E to an approximation of its integral (6). Realization of these operations is described here for completeness. ( i ) Let ( ,..., ) TN g g g be the vector of function values at the mesh nodes. If ( ) ( ) N i ii g y T y is the interpolating expansion in Chebyshev polynomials of the first kind, then A γ g , where the i -th column of the N N matrix A is the vector ( ( ),..., ( )) Ti i N
T y T y . Efficient calculation of this matrix is based on the recurrent relation: T T y , and k k k T yT T for k . Multiplication by this matrix realizes transition from the expansion coefficients to function values; the inverse matrix is used for transition from function values to the expansion coefficients. Similarly, interpolation by the expansion in Chebyshev polynomials of the second kind, ( ) ( ) N i ii g y U y , is related to the matrix B with the i -th column ( ( ),..., ( )) Ti i N
U y U y ; the matrix calculation is based on the relations:
1, 2
U U y , and k k k U yU U for k . ( ii ) For a given vector E of ( , ) k k E E t y values, the derivatives , / | k k t y E E y , see (14), are approximated by derivatives of the interpolating expansion ( ) N i ii
E T y with A γ E . Since d / d k k T y kU , we have ( 1) ( ) N i ii
E i U y . Let C be the N N matrix with elements
10 1 ij i j iC j i . Then BC D
E γ E , where the
N N differentiation matrix D BCA . ( iii ) Denote, for i N i i iT y y i i For ( ) N i ii
E T y we obtain d ( ) ( ) T T T T T T
E y A A γ σ E σ E σ E η , where ( ,..., ) TN σ and ( ) T T A η σ . We note that matrices , , , , A A B B D are used on each time step but, as well as the vector η , should be calculated only once. For, e.g., N this takes about 0.01 second and ensures the efficiency of our simulations. Appendix B. Computing
PM, / z H : interpolation from a lookup table We define a uniform mesh with ( 1) N nodes i in [0, 2 ] and find PM H analytically at the strip grid points ,..., N y y for every PM rotation angle i using (1). The lookup tables of PM,y ( , ) i k
H y and PM, ( , ) z k
H y are needed if the critical current density depends on magnetic field. The normal to strip component PM, ( , ) z k
H y is, for each k y , a periodic function of sampled at the angles i and PM, ( , ) / z i k
H y values were obtained by numerical differentiation in the Fourier space using the fast Fourier transform. In our simulations we chose N ; for N computing these tables takes about 0.1 second. Linear interpolation from the lookup tables provides a fast and accurate approximation to PM, ( , ) y k
H y , PM, ( , ) z k
H y , and PM, ( , ) / z k
H y for any rotation angle ( ) t ..