A new benchmark problem for electromagnetic modelling of superconductors: the high-Tc superconducting dynamo
Mark Ainslie, Francesco Grilli, Loic Queval, Enric Pardo, Fernando Perez-Mendez, Ratu Mataira, Antonio Morandi, Asef Ghabeli, Chris Bumby, Roberto Brambilla
AA new benchmark for electromagnetic modelling of superconductors: the high ‐ T c superconducting dynamo Mark
Ainslie , Francesco
Grilli , Loïc
Quéval , Enric
Pardo , Fernando
Perez ‐ Mendez , Ratu
Mataira , Antonio
Morandi , Asef
Ghabeli , Chris
Bumby , Roberto
Brambilla Department of Engineering,
University of Cambridge,
United
Kingdom Institute for
Technical
Physics,
Karlsruhe
Institute of Technology,
Karlsruhe,
Germany Group of Electrical
Engineering
Paris (GeePs),
CentraleSupélec,
University of Paris ‐ Saclay,
France Institute of Electrical
Engineering,
Slovak
Academy of Sciences,
Slovakia Robinson
Research
Institute,
Victoria
University of Wellington,
New
Zealand University of Bologna,
Bologna,
Italy Retired;
Ricerca sul
Sistema
Elettrico,
Milano,
Italy (formerly)
Corresponding author: [email protected]
Abstract
The high ‐ T c superconducting (HTS) dynamo is a promising device that can inject large DC supercurrents into a closed superconducting circuit. This is particularly attractive to energise HTS coils in NMR/MRI magnets and superconducting rotating machines without the need for connection to a power supply via current leads. It is only very recently that quantitatively accurate, predictive models have been developed which are capable of analysing these devices and explain their underlying physical mechanism. In this work, we propose to use the HTS dynamo as a new benchmark for the HTS modelling community.
The benchmark geometry consists of a permanent magnet (PM) rotating past a stationary HTS coated ‐ conductor wire (hereafter, HTS wire) in the open ‐ circuit configuration, assuming for simplicity the (infinitely long) case. Despite this geometric simplicity the solution is complex, comprising time ‐ varying spatially ‐ inhomogeneous currents and fields throughout the superconducting volume. In this work, this benchmark has been implemented using several different methods, including H ‐ formulation ‐ based methods, coupled H ‐ A and T ‐ A formulations, the Minimum
Electromagnetic
Entropy
Production method, and integral equation and volume integral equation ‐ based equivalent circuit methods. Each of these approaches show excellent qualitative and quantitative agreement for the open ‐ circuit equivalent instantaneous voltage and the cumulative time ‐ averaged equivalent voltage, as well as the current density and electric field distributions within the HTS wire at key positions during the magnet transit. Finally, a critical analysis and comparison of each of the modelling frameworks is presented, based on the following key metrics: number of mesh elements in the HTS wire, total number of mesh elements in the model, number of degrees of freedom (DOFs), tolerance settings and the approximate time taken per cycle for each model. . Introduction
The high ‐ T c superconducting (HTS) dynamo [1 ‐ is a promising device that can inject large DC supercurrents into a closed superconducting circuit. It could be used, for example, to energise HTS coils in NMR/MRI magnets and superconducting rotating machines without the need for connection to a power supply via current leads [4, Despite the extensive experimental work carried out to date, comprehensively understanding the underlying physical mechanism of such dynamo ‐ type flux pumps has proved challenging. A number of different explanations have been proposed to explain this mechanism [6 ‐ but quantitatively accurate, predictive calculations have been difficult. It was shown recently in Mataira et al . [14] that the open ‐ circuit voltage can be explained well – most importantly, with good quantitative agreement – using classical electromagnetic theory. The DC output voltage obtained from an HTS dynamo arises naturally from a local rectification effect caused by overcritical eddy currents flowing within the HTS wire [6 ‐ a classical effect that has been observed in HTS materials as far back as Vysotsky et al . [16]. The gap dependence of the open ‐ circuit voltage computed by Ghabeli and
Pardo [17] also agrees with experiments. In [17], it is also shown that this voltage is independent of the critical current density, J c , when the superconductor is fully penetrated by supercurrents. Since these overcritical eddy currents must recirculate within the
HTS wire, and can co ‐ exist with a transport current, the wire width is a key parameter and [18] shows that this should be sufficiently large so that the eddy and transport currents do not drive the full width of the stator into the flux ‐ flow regime. A number of different numerical models have now been developed as useful and cost ‐ efficient tools to further examine and explain the experimental results, as well as optimise and improve flux pump design. To adequately compare these, we propose a new benchmark for the HTS modelling community [19]: the
HTS dynamo. In this work, this benchmark problem is implemented using several different methods: • Coupled H ‐ A formulation [20]; • H ‐ formulation + shell current [14, • Segregated H ‐ formulation [22]; • Minimum
Electromagnetic
Entropy
Production (MEMEP) [23, • Coupled T ‐ A formulation [25, • Integral equation [27]; • Volume integral equation ‐ based equivalent circuit [28]. Section details the benchmark problem, including the geometry, parameters and the relevant assumptions. Section describes each of the modelling frameworks, including how the open ‐ circuit voltage is defined. Section presents the open ‐ circuit equivalent instantaneous voltage waveforms and the cumulative time ‐ averaged equivalent voltages, as well as the current density and electric field distributions within the HTS wire during the magnet transit.
Finally, in section a critical analysis and comparison of each of the modelling frameworks is presented, based on the following key metrics: number of mesh elements in the HTS wire, total number of mesh elements in the model, number of degrees of freedom (DOFs), tolerance settings and the approximate time taken per cycle for each model. . The
HTS dynamo benchmark
The geometry of the HTS dynamo benchmark problem is shown in Figure assuming for simplicity the (infinitely long) case. The permanent magnet (PM) rotates anticlockwise past the stationary
HTS wire in the open ‐ circuit configuration [14]. The PM has a width a and height b and a remanent flux density B r . The initial position of the PM is such that the centre of its outer face is located at (0, – R rotor ), i.e., θ M ( t = = – π /2. The
HTS wire has a width e and thickness f and is positioned such that its inner face is located (0, R rotor + airgap ). Table I lists the assumed parameters for the model, which are based on the model presented in [14] and correspond to the experimental setup in [29]. For simplicity, only the superconducting layer of the HTS wire is modelled and J c is assumed to be constant (where J c = I c /( e ꞏ f )). It was shown in [14] that this assumption does not impact the essential dynamics to deliver a DC voltage, which is simply that the wire must exhibit a non ‐ linear resistivity. Isothermal conditions are assumed (i.e., a constant temperature, T ) and hence no thermal model needs to be included. Regardless of the modelling framework used, the open ‐ circuit equivalent instantaneous voltage and the cumulative time ‐ averaged equivalent voltage waveforms shown later in Figures and respectively, should be obtained by implementing the benchmark model. Figure Geometry of the HTS dynamo benchmark. A permanent magnet (PM) rotates anticlockwise past an HTS wire. able I. HTS dynamo benchmark parameters.
Permanent magnet (PM)
Width, a mm Height, b mm Active length (depth), L mm Remanent flux density, B r T HTS stator wire
Width, e mm Thickness, f μ m Critical current, I c [self ‐ field, K] A n value Rotor external radius, R rotor mm Distance between PM face & HTS surface, airgap mm Frequency of rotation Hz Number of cycles Modelling frameworks
General definitions
The nonlinear resistivity, ρ ( J ), of the superconductor is simulated using an E ‐ J power law [30 ‐ nc c E JJ J E J (1) where, in Cartesian coordinates and infinitely long (2D) problems in the z direction, J = [0 J z ] and E = [0 E z ] are the current density and electric field, respectively, which are assumed to be parallel to each other such that E = ρ J . E = μ V/cm is the characteristic electric field and n defines the steepness of the transition between the superconducting state and the normal state. For n > equation (1) becomes a reasonable approximation of the Critical
State
Model (CSM), for which n approaches infinity [33, although accurate agreement with the CSM may require n values in the range of ‐ [24]. In general, the instantaneous measured voltage, V ( t ), is the path integral of the gradient of the electrostatic potential, ∇ φ , over the superconductor and measurement wires [14, When the excitations are periodic with period T (external magnetic field and transport currents) and for infinitely long geometries, the DC component of the voltage DC t Tt V V t dtT (2) corresponds to [14, C ( ') ' t T avet LV E t dtT (3) where L is the active length of the dynamo, i.e., the active length (depth) of the PM, and E ave is the electric field, E z , averaged over the cross ‐ section of the superconductor, S :
1( ) ( , , ) ave zS
E t E x y t dSS (4) Then, the equivalent instantaneous voltage, V eq ( t ), is defined as ( ) ( ) eq ave V t L E t (5) and the cumulative time ‐ averaged equivalent voltage, V cumul ( t ), is given by
1( ) ( ) tcumul eq
V t V t dtt (6) which corresponds to V DC for large enough t. Under open ‐ circuit conditions, no net transport current flows, such that, at all times ( ) 0 zS I t J dS (7) which is implemented as a constraint in each of the models. H ‐ formulation models For the H ‐ formulation [36 ‐ the independent variables are the components of the magnetic field strength, H = [ H x H y and the governing equations are derived from Maxwell’s equations – namely, Ampere’s (8) and
Faraday’s (9) laws:
H J (8) t BE = (9)
The permeability µ = µ , and equations (8) and (9) are combined with the E ‐ J power law, equation (1). Coupled H ‐ A formulation (H ‐ A) The coupled H ‐ A formulation, proposed for modelling superconducting rotating machines in [20], models the entire rotating model, with the H ‐ formulation solved in a small region local to the HTS wire and the magnetic vector potential A solved elsewhere (thus, much of the model follows the usual construction for conventional rotating machines). In such a mixed ‐ formulation model, careful attention must be paid to coupling variables across common boundaries between the H and A ubdomains to maintain continuity: this is achieved by coupling, in weak form, the electric field from the A ‐ formulation to the H ‐ formulation and coupling the tangential components of the magnetic field from the H ‐ formulation to the A ‐ formulation, equivalent to a Neumann boundary condition [20]. In the implementation of this model here, a simplification is made by limiting the region that directly solves the vector fields associated with Maxwell’s equations to a small region surrounding the conductive (current ‐ carrying) subdomain, i.e., the H ‐ formulation subdomain including the HTS wire.
This allows most of the model to use the magnetic scalar potential, V m , to calculate the PM field, for which the following magnetic flux conservation equation holds: m r V B (10) where B r is the remanent flux density (1.25 T for the PM assumed here – see Table I – and zero elsewhere). In the magnetic vector potential formulation, the magnetic flux density is defined as B A (11) and the electric field as t AE (12) automatically fulfilling Faraday’s law (equation (9)) and the magnetic flux conservation law B (13) and then Ampere’s law (equation (8)) is solved. This model is implemented in COMSOL
Multiphysics® using the
Rotating
Machinery,
Magnetic (RMM) interface in the AC/DC module, which uses this mixed formulation of V m and A . An appropriate gauge is chosen such that the scalar electric potential (see equations (22) and (37)) vanishes and only A has to be considered (equation (12)). The H ‐ formulation is implemented in the Magnetic
Field
Formulation (MFH) interface, also in the AC/DC module. In addition to the H ‐ A coupling described above, the in ‐ built Mixed
Formulation
Boundary node in the RMM interface imposes continuity between V m and A on either side of the magnetic scalar/vector potential boundary, such that m V m V A n H n H n (14) and m V A n B n B n A (15) where the surface normals, n and n , are antiparallel ( n = – n ). Equations (14) and (15) are implemented as weak contributions and can be interpreted as a surface current density, s J n H , and magnetic surface charge density, m = n B , respectively. .2.2. H ‐ formulation + shell current (H+SC) Modelling the rotation of the PM in Fig. is not trivial in the H ‐ formulation. The approach developed in [14], and used again in [18, is to represent the PM as a time ‐ dependent sheet current K sheet along the boundary of the rotor domain ∂ Ω R . This reproduces the rotating field of the PM in the domain of interest (i.e., outside the rotor), while not introducing a rotating mesh or inter ‐ model couplings to capture the rotation. The major advantage of this method is that it does not require the self ‐ field correction described in Section (segregated H ‐ formulation), and ensures the whole model is solved natively as a finite element problem. However, the major disadvantage of this is the large number of mesh elements committed to simulating the boundary of the rotor domain which leads to a higher computational cost. To implement the shell model, the PM is simulated in the sub ‐ domain of the rotor Ω R in a static model. Setting the boundary condition of the model to be magnetically insulating along the boundary of the rotor ∂ Ω R : R Ω n E (16) This gives a solution where the flux of the magnet is completely enclosed in the rotor domain and H = elsewhere. By the principle of superposition, it must be the case that the field of the PM is contained inside the boundary ∂ Ω R by a shell current, K shell , on this boundary, that produces the opposite PM field outside the boundary, i.e., H m + H shell = Hence, the effect of the shell current is to produce the image of the PM’s magnetic field outside the boundary of the initial model. The full problem can now be solved by omitting the original magnet and instead using the negative of the shell current distribution: ( ( )) R shell M Ω t H K (17) where θ is the angular coordinate around the rotational axis of the rotor. This produces the magnetic field of the PM, in the domain of interest, and rotates it by the angle θ M ( t ) with time. Segregated H ‐ formulation (SEG ‐ H) The segregated model is comprised of a magnetostatic PM model and a time ‐ dependent H ‐ formulation HTS wire model.
The former is coupled unidirectionally to the latter using boundary conditions [22] and a translation (rotation) operator for the PM’s static magnetic field (see
Figure
This avoids, like the preceding shell current model, the need to model the rotating PM (e.g., using a moving mesh), but also significantly reduces the number of mesh elements in the HTS model.
Figure Segregated finite ‐ element model: (a) magnetostatic PM model, (b) time ‐ dependent H ‐ formulation HTS wire model. On the outer boundary of the H ‐ formulation model, the sum of the applied field, H ext , and the self ‐ field, H self , is applied as a Dirichlet boundary condition. To mimic the rotation, H ext is obtained by rotating the field of a static permanent magnet, H PM : ext,x PM,x rot rotM Mext,y PM,y rot rotM M ( , , ) ( , )cos ( ) sin ( )( , , ) ( , )sin ( ) cos ( ) H x y t H x yt tH x y t H x yt t (18) rotrot cos ( ) sin ( )sin ( ) cos ( )
M MM M x t t xy t t y (19) where θ M ( t ) is the rotor angle and ( x rot , y rot ) are the coordinates in the rotated coordinate system. The self ‐ field, H self , created by the supercurrent flowing in the HTS wire, is calculated at each time step by numerical integration of the Biot ‐ Savart law over the
HTS wire subdomain: zself,x 2 2 ( ', ', ) ( ')1( , , ) ' '2 ( ') ( ') S J x y t y yH x y t dx dyx x y y (20) zself,y 2 2 ( ', ', ) ( ')1( , , ) ' '2 ( ') ( ') S J x y t x xH x y t dx dyx x y y (21) Minimum
Electromagnetic
Entropy
Production (MEMEP)
MEMEP is a variational method that has been shown to be ideally suited modelling materials with highly nonlinear E ( J ) relationships, such as superconductors [23, MEMEP solves the current density J by minimizing a functional that contains all the variables of the problem, including the magnetic vector potential A , current density J and scalar potential φ . It has been proven that this functional always presents a minimum, it is unique, and it is the solution of Maxwell’s equations in differential form [24]. his method is fast because the current density only exists inside the superconducting region, and thus the mesh is only required inside this region. The general equation for the current density and the scalar potential are t AE (22) J (23) In Coulomb’s gauge ( A ), A can be separated into the contributions from the current density, A J , and the external sources, A a , [23]. For infinitely long problems in the z direction (or J = [0 J z ], E = [0 E z ], and A = [0 A z ]. The A J contribution in Coulomb’s gauge follows [33] ( ) ( ) ln2 J S
A r dS'J r' r r' (24) Equation (23) is always satisfied for problems, and thus only equation (22) needs to be solved. To solve this equation, the following functional should be minimised [23, J aS
A AL ds J J U J Jt t (25) where U is the dissipation factor defined as [24] ( ) ( ) J U J E J' dJ' (26) This dissipation factor can include any E ‐ J relationship for superconductors, including the multi ‐ valued relation of the CSM [24, In this problem, the non ‐ uniform applied magnetic field caused by the rotating PM appears in the functional in the form of A a . Then, the program only requires calculation of the vector potential once for each time step within the first cycle. The impact on the total computing time is negligible because the minimisation takes most of the computing time. The vector potential generated by the PM with uniform magnetisation M can be calculated by the magnetization sheet current density K = M × e n , where M is the PM magnetisation and e n is the unit vector normal to the surface. For uniform magnetisation, the vector potential A M generated by the PM is: ( ) ( ) ln2 M m nS r M dl r r r' A e e (27) where e m is the unit vector in the magnetisation direction, ∂ S represents the edges of the PM cross ‐ section, and dl ’ is the length differential on the edge. In this work, A M was evaluated numerically. Note that the cross product in the equation above always follows the z direction, since both e m and e n are in the xy plane. MEMEP can also take a J c ( B , θ ) dependence into account by solving J , then calculating B , and iterating until the difference is below a certain tolerance [23]. In this case, B from the PM should also be calculated, as carried out in [17]. .4 Coupled T ‐ A formulation (T ‐ A) The T ‐ A formulation was proposed in [25, to tackle the problem of simulating superconductors characterised by a very high aspect (width:thickness) ratio, such as HTS wires.
The main idea is to use the magnetic vector potential A for calculating the magnetic field in the whole domain and the current vector potential T for calculating the current density J T (28) in the superconductor. The obtained current density is re ‐ injected as an external current density in the A ‐ formulation part. The superconductor can be simulated as a object in this problem, the line representing the superconducting wire. Further resulting simplifications are that the current vector potential has only one component and that the transport current flowing in the superconductor can be imposed by means of simple boundary conditions for T at the wire’s edges. In particular, the current is determined by the difference between the values of T at the wire’s edges. In the benchmark considered here, there is no transport current and the simple condition T = any constant is applied at each of the wire’s edges. The T ‐ A formulation has been recently extended to simulate superconductors of finite thickness [46]. In this case, T has two components and the setting of the boundary conditions is a little less intuitive – see [46] for details. To distinguish between these two T ‐ A formulations, we refer to these hereafter as T ‐ A (1D) and T ‐ A (2D), respectively. In addition, these two formulations are implemented separately using only the magnetic vector potential A (implemented in COMSOL using the
Magnetic
Field (MF) interface), referred to as VP (vector potential), and using the mixed scalar ‐ vector potential detailed in Section implemented using the
Rotating
Machinery,
Magnetic (RMM) interface and referred to as SP (scalar potential). Integral equation (IE)
The current distribution along a segment representative of an infinitely thin (1D) superconducting layer (as per Section can be given by an integral equation that can be easily solved by the finite ‐ element method. This approach was proposed in [27], and later extended to consider interacting HTS wires in [47, For the problem analysed here, the integral equation is written in the Partial
Differential
Equation (PDE) module of COMSOL in and takes the following form ( ) s J f Q K C (29) with ( , ) ( , ) x t na
K x t H t d (30) and
1( , ) ( , ) ln2 a t sa
Q x t J t x d (31) where ρ is the power ‐ law resistivity (equation (1)), J s is the sheet current density (current per unit width, A/m) in the z direction, f is the thickness of the superconductor (see Figure a is the half ‐ idth (i.e., e = a ) and H n is the normal component of the external field impinging on the superconductor. The constant C is set at each time step to satisfy the constraint on the desired transport current, which in this benchmark is always zero (see equation (7)): ( , ) 0 a sa J t d (32) To obtain the external field, the PDE module is coupled to the Magnetic
Field (MF) interface in the AC/DC module, which calculates only the magnetic field generated by the rotating PM.
The mesh rotation is considered by using the arbitrary Lagrangian ‐ Eulerian (ALE) method.
The MF module does not contain the “reaction” term of the field created by the currents flowing in the superconductor. As a consequence, in order to visualize the total magnetic field in the whole simulated domain, one needs to use a second MF module where the magnetic field generated by the sheet current is added to the one generated by the permanent magnet. This has, of course, an additional computational cost because of the additional DOFs.
For the benchmark proposed here, where the focus is the quantities in the superconductor, this second MF module was not used. The main advantage of this method relies in the fact that the current density is the state variable of the equation and, as a consequence, it is not obtained by the spatial derivative as required by Ampere’s law (equation (8)).
Note that in equations (29) ‐ (31), the time derivative of J S appears under the integral sign, but using standard procedure (Carleman’s equation), it is possible reformulate these in a way that extracts the time derivative: ( , )2( , ) ( ) at s a F t aJ x t d C tf x a x (33) where ( , ) ( , ) ( , ) s F x t J x t fK x t (34)
Therefore, ( , )2( , ) ( ) at s a
F t aJ x t d C tx a x (35) where ( , )( , ) ( ) ( , ) ( , ) ( , ) s E tF t J t K t K tf (36)
This second form may be more appropriate for standard numerical routines for solving integral equations (Nystroem method) that do not require finite elements on the strip segment. Volume integral equation ‐ based equivalent circuit (VIE) The volume integral equation ‐ based equivalent circuit is obtained by separating the total electromotive force at any point of the superconductor into two contributions: one contribution due to the time ‐ varying field produced by the current induced in the superconductor and a second contribution due to the movement of the PM.
According to this, and by expressing the magnetic flux density via the magnetic vector potential as in equation (11), Faraday’s law at any point in the superconductor gives [49 ‐ nt PM t AE v B (37) where A int is the vector potential of the current in the superconductor, B PM is the PM field and v is the velocity of B PM at the considered point, expressed in the fixed reference frame of Fig. A numerical solution of the problem is obtained by subdividing the superconductor domain into a finite number of elements and by enforcing equation (37) to be satisfied, in the weak form, over each element of the discretization [28, The state variables of the problem are the current densities of each element. This is obtained by relating E and A int in equation (37) to the current density via the E ‐ J power law (equation (1)) and equation (24), respectively. The thin shell model (see
Sections and is used for calculating the field of the PM at each position. Each of the discretised equations obtained via the weighted residual approach corresponds to the voltage balance of a circuit branch involving a non ‐ linear resistor arising from the electric field in the superconductor, a coupled inductor representing the magnetic interaction of the induced current, and a voltage generator corresponding to the Lorentz ‐ like electromotive term. Hence, this circuit picture gives rise to the name of the method. Results
Figure shows a comparison of the open ‐ circuit equivalent instantaneous voltage waveforms calculated by each of the models for the nd transit of the PM past the HTS wire, ignoring any initial transient effects that may be present in the st cycle. Indeed, qualitatively, the distinct four peaks and noticeable left ‐ to ‐ right asymmetry observed in experiments (see [14], for example) are reproduced and there is excellent quantitative agreement between the models for the magnitude of these peaks. The particular characteristics of this voltage waveform give rise to the DC output voltage of the HTS dynamo, which can be further evidenced by examining the cumulative time ‐ averaged equivalent voltage, V cumul , given by equation (6). V cumul calculated by each of the models over the cycles of PM rotation is shown in Figure In all cases this converges to a non ‐ zero asymptotic value, clearly showing the DC output. Again there is excellent qualitative and quantitative agreement between each of the models: the average value of V cumul is –9.41 µV with a standard deviation of µV. It should be noted that, as described earlier, there is some discrepancy between these results and those observed in experiments because of the use of the constant J c approximation. The use of the angular field ‐ dependence of J c ( B , θ ) to consider the suppression of J c with magnetic field, as considered elsewhere in [14, is needed for good agreement with experiment and can be done easily by modifying the E ‐ J power law, equation (1), such that J c = J c ( B , θ ). Figure shows the current density normalised to J c0 , J / J c0 , and electric field, E , distributions within the wire for three key PM positions as the magnet travels past it: as the magnet approaches the centre of the wire from the right ‐ hand side ( t = ms in the nd cycle); when the magnet is aligned with the centre of the wire ( t = ms); and as the magnet moves away from the wire on its left ‐ hand side ( t = ms). The dynamics of the current flowing within the wire and the related local electric fields ultimately give rise to the voltage waveforms shown in Figure with close to the first negative peak, close o the trough in between the nd and rd (positive) peaks and close to the fourth negative peak. At an overcritical ( J > J c ) eddy current flows in the right ‐ hand side of the wire, which then returns (i.e., flows in the opposite direction) on the other side of the wire at a lower magnitude, giving rise to the asymmetric electric field distribution as described by the E ‐ J power law. At these forward and reverse currents are almost equal and opposite, such that V ≈ and at the reverse situation to occurs, such that the J and E profiles are essentially mirrored. Similar profiles were also obtained in [14] for the J c ( B , θ ) case, except that the suppression of J c with magnetic field results in higher local electric fields in the region close to the PM.
Figure Open ‐ circuit equivalent instantaneous voltage, V eq ( t ), waveforms calculated by each of the models for the nd transit of the PM past the HTS wire, ignoring any initial transient effects that may be present in the st cycle. Qualitatively, the distinct four peaks and noticeable left ‐ to ‐ right asymmetry observed in experiments are reproduced and there is excellent quantitative agreement between the models for the magnitude of these peaks. Figure Cumulative time ‐ averaged equivalent voltage, V cumul ( t ), calculated for each of the models over the cycles of PM rotation, clearly showing a DC output, with excellent qualitative and quantitative agreement between each of the models. The average value of V cumul after cycles of PM rotation is –9.41 µV with a standard deviation of µV. Figure Distributions within the wire for: current density normalised to J c0 , J / J c0 , and electric field, E . The calculated distributions are shown for three key PM positions: as the magnet approaches the centre of the wire from the right ‐ hand side ( t = ms in the nd cycle); when the magnet is aligned with the centre of the wire ( t = ms); and as the magnet moves away from the wire on its left ‐ hand side ( t = ms). The dynamics of the current flowing within the wire and the related local electric fields ultimately give rise to the voltage waveforms shown in Figure . Discussion As shown in Figures ‐ all of the models produced the expected benchmark solution with excellent qualitative and quantitative agreement. In this section, a critical analysis and comparison of each of the modelling frameworks is presented. Table II lists the key metrics assessed for each benchmark model: the number of mesh elements in the HTS wire; the total number of mesh elements in the model; the number of DOFs; the relative and absolute tolerance settings (for
FEM ‐ based models), tolerance for the mutual interaction matrix (MEMEP, programmed in C++) or ode23b solver relative tolerance (VIE, programmed in MATLAB); and the approximate time taken per cycle for each model. In the interest of a fair comparison, all of the FEM ‐ based models were run on the same computer under the same conditions (e.g., COMSOL
Multiphysics™ version and, where possible, the number of mesh elements in the HTS wire set to x along the width and thickness, respectively. The key findings and comparisons are detailed below: The clear winner in terms of computational speed is the MEMEP method, with the entire cycles taking a little over two minutes to solve. This can be explained by the limited number of DOFs because only the
HTS wire needs to be meshed. It should also be noted that the MEMEP model was run on a slightly inferior processor, so slightly overestimates the computational time per cycle. The next best performers are the
SEG ‐ H and VIE methods, which are also modelling frameworks that emphasise a reduced number of mesh elements. It should be noted that, in the FEM ‐ based models, the mesh was optimised as best possible as a compromise between accuracy and computational speed. Thus, there is scope in many of the FEM ‐ based models to improve their accuracy somewhat by using a finer mesh, but there will be an associated increase in computational time. Several models (H ‐ A, T ‐ A (2D), SEG ‐ H and H+SC) took advantage of the artificial expansion technique presented in [53, to increase the HTS layer thickness from µm to µm, improving the computational speed without compromising accuracy. The rotating machine ‐ like models that made use of the mixed scalar ‐ vector potential (H ‐ A, T ‐ A (1D) SP and T ‐ A (2D) SP) also performed well in terms of computational time, and the integral equation ‐ based model (IE), even with the use of the magnetic vector potential with second ‐ order (quadratic) elements and the associated significant increase in DOFs, also performed comparably. The T ‐ A (2D) formulations – both SP and VP – are found to be reasonably unstable, and even a finer mesh (60 x and tighter tolerances settings could not improve this performance dramatically. The spikes seen in the J and E plots in Figure are from the T ‐ A (2D) VP model; the E spikes are much more pronounced due to E being proportional to J n (see equation (1)). Spurious oscillations of a similar kind were also presented in [45, for the T ‐ A (1D) model used for AC loss calculations. It was recommended in both [45, that second ‐ order (quadratic) elements be used for A to mitigate this, which does have its associated computational cost. However, it should also be noted that for the benchmark here and in [45, this did not impact the calculations of interest significantly (AC loss and voltage, respectively), although the voltage waveforms in Figure are clearly noisier in comparison to other models. The use of the scalar potential (H ‐ A, T ‐ A (1D) SP and T ‐ A (2D) SP) – with first ‐ order (linear) elements – improved stability and spurious oscillations/noise, as well as computational speed, in comparison to the VP models. Using this mixed scalar ‐ vector formulation is a potentially useful alternative for modelling such dynamos, as well as superconducting rotating machines in general. For models in particular, the scalar potential formulation ntroduces fewer DOFs and can ensure a more accurate coupling of the magnetic field. Indeed, although all of these models were specifically created for the HTS dynamo benchmark, many of the findings are equally applicable to and useful for modelling superconducting rotating machines. In terms of ease of use, all models (except for the MEMEP and
VIE methods, which are self ‐ programmed using C++ and
MATLAB, respectively) were implemented in COMSOL
Multiphysics™.
COMSOL is a popular commercial software package with a reasonably gentle learning curve and is currently used by dozens of research groups worldwide to model superconductivity ‐ related problems [56, Many shared modelling examples, most associated with peer ‐ reviewed publications, are available on the HTS
Modelling
Workgroup website [58].
COMSOL now has a dedicated superconductivity interface (the MFH interface; see
Section with dedicated technical support. However, there is a significant disadvantage when using commercial software in that users do not have complete control over its implementation: some of the programming cannot be accessed easily, if at all. It is of particular note that the SEG ‐ H and H+SC models, which were built and optimised in COMSOL version ran significantly slower when opened and run in version The computational time per cycle for both versions is included in Table
II. At the time of writing, COMSOL were unable to explain why these models had longer run times in version despite backward compatibility. There is also an associated financial cost that can be a barrier to some researchers, which is where self ‐ programmed techniques and those implemented in free software have a distinct advantage. Table
II.
Key metrics assessed for each benchmark model.
Model
Mesh (SC)
Mesh (total)
DOFs
Rel./abs. tolerance
Approx. time/cycle [min/cycle]
Software implementation
MEMEP (120 x ‐ < A C++
SEG ‐ H (120 x ‐ / B COMSOL B COMSOL
VIE (120 x ‐ / ‐ B MATLAB H ‐ A (120 x ‐ / B COMSOL T ‐ A (2D) SP (60 x ‐ / ‐ B COMSOL IE (120 x ‐ / B COMSOL T ‐ A (1D) SP (120 x ‐ / ‐ B COMSOL
H+SC (120 x ‐ / ‐ B COMSOL > COMSOL T ‐ A (1D) VP (120 x ‐ / B COMSOL T ‐ A (2D) VP (60 x ‐ / B COMSOL PC specifications: A Intel®
Core™ i7 ‐ CPU @ GHz, GB RAM (10% memory used for
MEMEP model),
Ubuntu
LTS, ‐ bit B Intel®
Core™ i9 ‐ CPU @ GHz, GB RAM,
Microsoft
Windows Pro, ‐ bit Other notes: Tolerance for the mutual interaction matrix Default settings for
MATLAB/ ode23b solver . Conclusion In this work, a new benchmark problem for the HTS modelling community – the HTS dynamo – was proposed, consisting of a permanent magnet (PM) rotating past a stationary HTS wire in the open ‐ circuit configuration. The benchmark was then implemented using several different methods, including H ‐ formulation ‐ based methods, coupled H ‐ A and T ‐ A formulations, the Minimum
Electromagnetic
Entropy
Production method, and integral equation and volume integral equation ‐ based equivalent circuit methods. Excellent qualitative and quantitative agreement was obtained between all models for the open ‐ circuit equivalent instantaneous voltage and the cumulative time ‐ averaged equivalent voltage, as well as the current density and electric field distributions within the HTS wire at key positions during the magnet transit. The average value for all the models of the DC output voltage of the HTS dynamo, determined by the cumulative time ‐ averaged equivalent voltage over cycles of PM rotation, was calculated to be –9.41 µV with a standard deviation of µV. A critical analysis and comparison of each of the modelling frameworks was presented, based on the following key metrics: number of mesh elements in the HTS wire, total number of mesh elements in the model, number of DOFs, tolerance settings and the approximate time taken per cycle for each model.
The clear winner in terms of computational speed is the MEMEP method, with the entire cycles taking around two minutes to solve, due to the limited number of DOFs because only the
HTS wire needs to be meshed. The next best performers were the
SEG ‐ H and VIE methods, which are also modelling frameworks that emphasise a reduced number of mesh elements. Several models took advantage of an artificial expansion technique to increase the HTS layer thickness from µm to µm, improving the computational speed without compromising accuracy. A number of models use a rotating machine ‐ like modelling framework – in particular, the coupled H ‐ A and T ‐ A formulations – and it is shown that the use of a mixed scalar ‐ vector potential (implemented using COMSOL’s
Rotating
Machinery,
Magnetic interface) results in a significant improvement in both computational speed and stability, compared to models that use only the vector potential (implemented using COMSOL’s
Magnetic
Field interface). In the latter case, it is recommended to use second ‐ order (quadratic) elements for A – in particular, for the T ‐ A formulation – to mitigate against spurious oscillations and improve stability, which has an associated computational cost. Using the mixed scalar ‐ vector formulation provides a potentially useful alternative for modelling such dynamos, as well as superconducting rotating machines in general. This benchmark and the results contained herein now provide researchers with a suitable framework to validate, compare and optimise their own methods for modelling the HTS dynamo.
Acknowledgements MA would like to acknowledge financial support from an Engineering and
Physical
Sciences
Research
Council (EPSRC)
Early
Career
Fellowship
EP/P020313/1.
CB, RM and MA were supported in part by NZ Royal
Society
Marsden
Grant
No.
MFP ‐ VUW1806. FP would like to acknowledge the support of the Consejo
Nacional de Ciencia y Tecnologia and
Secretaria de Energia de Mexico,
Reference
No.
Additional data related to this publication are available at the University of Cambridge data repository (https://doi.org/10.17863/CAM.54005).
RCID iDs
Mark
Ainslie ‐ ‐ ‐ Francesco
Grilli ‐ ‐ ‐ Loïc
Quéval ‐ ‐ ‐ Enric
Pardo ‐ ‐ ‐ Fernando
Perez ‐ Mendez ‐ ‐ ‐ Ratu
Mataira ‐ ‐ ‐ Antonio
Morandi ‐ ‐ ‐ Asef
Ghabeli ‐ ‐ ‐ Chris
Bumby ‐ ‐ ‐ Roberto
Brambilla
References [1]
Hoffmann C, Pooke D and Caplin A D IEEE
Trans.
Appl.
Supercond. ‐ [2] Volger J and Admiraal P Phys.
Lett. ‐ [3] Beelen H V et al Physica ‐ [4] Bai Z et al Cryogenics ‐ [5] Bumby C W et al Supercond.
Sci.
Technol. [6] Geng J et al Appl.
Phys.
Lett. [7]
Geng J et al J. Phys. D [8] Bumby C W et al Appl.
Phys.
Lett. [9]
Campbell A M Supercond.
Sci.
Technol. [10] Wang W and Coombs T Phys.
Rev.
Appl. [11] Giaever I IEEE
Spectrum ‐ [12] Kaplunenko V, Moskvin S and Schmidt V Fiz.
Nizk.
Temp. [available at: https://fnte.ilt.kharkov.ua/fnt/pdf/11/11 ‐ ‐ [13] van de Klundert L and ten Kate H Cryogenics [14] Mataira R C, Ainslie M D, Badcock R A and Bumby C W Appl.
Phys.
Lett. [15]
Jiang Z et al Appl.
Phys.
Lett. [16]
Vysotsky V S et al Supercond.
Sci.
Technol. ‐ [17] Ghabeli A and Pardo E Supercond.
Sci.
Technol. [18] Mataira R, Ainslie M D, Badcock R and Bumby C W IEEE
Trans.
Appl.
Supercond. [19] HTS
Modelling
Workgroup:
Benchmarks [20]
Brambilla R et al IEEE
Trans.
Appl.
Supercond. [21] Mataira R et al Phys.
Rev.
Appl. at press [22] Quéval L et al Supercond.
Sci.
Technol. [23] Pardo E, Ŝ ouc J and Frolek L Supercond.
Sci.
Technol. [24] Pardo E and Kapolka M J. Comput.
Phys. ‐ [25] Zhang H et al Supercond.
Sci.
Technol. [26] Benkel T et al IEEE
Trans.
Appl.
Supercond. [27] Brambilla R et al Supercond.
Sci.
Technol. [28] Morandi A and Fabbri M Supercond.
Sci.
Technol. [29] Badcock R A et al IEEE
Trans.
Appl.
Supercond. [30] Plummer C J G and Evetts J E IEEE
Trans.
Magn. ‐ [31] Rhyner J Physica C ‐ [32] Brandt E H Phys.
Rev. B ‐ [33] Ainslie M D and Fujishiro H IOP
Expanding
Physics (Bristol:
IOP
Publishing) [DOI: ‐ ‐ ‐ ‐ p ‐ [34] Grilli F et al IEEE
Trans.
Appl.
Supercond. [35] Clem J R Phys.
Rev. B ‐ [36] Kajikawa K et al IEEE
Trans.
Appl.
Supercond. ‐ [37] Pecher R et al Proc.
EUCAS pp [38] Hong Z, Campbell A M and Coombs T A Supercond.
Sci.
Technol. ‐ [39] Brambilla R, Grilli F and Martini L Supercond.
Sci.
Technol. ‐ [40] Ainslie M D, Flack T J, Hong Z and Coombs T A Int. J. Comput.
Math.
Electr.
Electron.
Eng. ‐ [41] Ainslie M D et al IEEE
Trans.
Appl.
Supercond. ‐ [42] Ainslie M D, Flack T J and Campbell A M Physica C ‐ [43] Bossavit A IEEE
Trans.
Magn. ‐ [44] Prigozhin L IEEE
Trans.
Appl.
Supercond. ‐ [45] Liang F et al J. Appl.
Phys. [46]
Grilli F et al [https://arxiv.org/abs/2004.01913] [47] Grilli F et al IEEE
Trans.
Appl.
Supercond. ‐ [48] Brambilla R et al Supercond.
Sci.
Technol. [49] Morandi A et al IEEE
Trans.
Appl.
Supercond. [50] Fabbri M et al IEEE
Trans.
Magn. ‐ [51] Perini E, Giunchi G, Geri M and Morandi A IEEE
Trans.
Appl.
Supercond. ‐ [52] Morandi A Supercond.
Sci.
Technol. [53] Hong Z and Coombs T A J. Supercond.
Nov.
Magn. ‐ [54] Zermeno V M R et al J. Appl.
Phys. [55]
Berrospe ‐ Juarez E, Zermeno V M R, Trillaud F and Grilli F Supercond.
Sci.
Technol. [56] Shen B, Grilli F and Coombs T Supercond.
Sci.
Technol. [57] Shen B, Grilli F and Coombs T IEEE
Access ‐ [58] HTS
Modelling
Workgroup:
Shared models