A new benchmark for electromagnetic modelling of superconductors: the high-Tc superconducting dynamo

The high-Tc 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 rotating past a stationary HTS coated-conductor wire in the open-circuit configuration, assuming for simplicity the 2D (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][2][3] 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,5]. 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][7][8][9][10][11][12][13], 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-8, 14, 15]: 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,18,21]; • Segregated H-formulation [22]; • Minimum Electromagnetic Entropy Production (MEMEP) [23,24]; • Coupled T-A formulation [25,26]; • Integral equation [27]; • Volume integral equation-based equivalent circuit [28].
Section 2 details the benchmark problem, including the geometry, parameters and the relevant assumptions. Section 3 describes each of the modelling frameworks, including how the open-circuit voltage is defined. Section 4 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 5, 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 1, assuming for simplicity the 2D (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 = 0) = -π/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 timeaveraged equivalent voltage waveforms shown later in Figures 3 and 4, respectively, should be obtained by implementing the benchmark model.  The nonlinear resistivity, ρ(J), of the superconductor is simulated using an E-J power law [30][31][32]: (1) where, in Cartesian coordinates and infinitely long (2D) problems in the z direction, J = [0 0 J z ] and E = [0 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 0 = 1 μ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 > 20, equation (1) becomes a reasonable approximation of the Critical State Model (CSM), for which n approaches infinity [33,34], although accurate agreement with the CSM may require n values in the range of 100-1000 [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,35]. 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 corresponds to [14,17,21] 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: Then, the equivalent instantaneous voltage, V eq (t), is defined as (5) and the cumulative time-averaged equivalent voltage, V cumul (t), is given by which corresponds to V DC for large enough t. Under open-circuit conditions, no net transport current flows, such that, at all times which is implemented as a constraint in each of the models.

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 subdomains 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: 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    A E (12) automatically fulfilling Faraday's law (equation (9)) and the magnetic flux conservation law 0   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 and (15) where the surface normals, n 1 and n 2 , are antiparallel (n 1 = -n 2 ). 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.

H-formulation + shell current (H+SC)
Modelling the rotation of the PM in Fig. 1 is not trivial in the H-formulation. The approach developed in [14], and used again in [18,21], 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 3.2.3 (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 : (16) This gives a solution where the flux of the magnet is completely enclosed in the rotor domain and H = 0 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 = 0. 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: (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 Hformulation 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 2). 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. On the outer boundary of the H-formulation model, the sum of the applied field, H ext , and the selffield, 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 : 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 2D Biot-Savart law over the HTS wire subdomain:

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,24]. 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].
This 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 In Coulomb's gauge ( , 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 2D), Equation (23) is always satisfied for 2D problems, and thus only equation (22) needs to be solved. To solve this equation, the following functional should be minimised [23,24]: where U is the dissipation factor defined as [24] This dissipation factor can include any E-J relationship for superconductors, including the multivalued relation of the CSM [24,43,44]. 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: (27) where e m is the unit vector in the magnetisation direction, ∂S represents the edges of the PM crosssection, 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].

Coupled T-A formulation (T-A)
The T-A formulation was proposed in [25,45] 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 1D object in this 2D problem, the 1D 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 scalarvector potential detailed in Section 3.2.1, 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 3.4) can be given by an integral equation that can be easily solved by the finiteelement method. This approach was proposed in [27], and later extended to consider interacting HTS wires in [47,48].
For the problem analysed here, the integral equation is written in the Partial Differential Equation (PDE) module of COMSOL in 1D and takes the following form 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 1), a is the half-width (i.e., e = 2a) 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)): To obtain the external field, the 1D PDE module is coupled to the 2D 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: Therefore, 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][50][51] (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. 1. A numerical solution of the problem is obtained by subdividing the superconductor domain into a finite number of 2D elements and by enforcing equation (37) to be satisfied, in the weak form, over each element of the discretization [28,52]. 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 3.2.2 and 3.3) 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. Figure 3 shows a comparison of the open-circuit equivalent instantaneous voltage waveforms calculated by each of the models for the 2 nd transit of the PM past the HTS wire, ignoring any initial transient effects that may be present in the 1 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 timeaveraged equivalent voltage, V cumul , given by equation (6). V cumul calculated by each of the models over the 10 cycles of PM rotation is shown in Figure 4. In all cases this converges to a non-zero asymptotic value, clearly showing the DC output.

Results
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 0.34 µ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 fielddependence of J c (B, θ) to consider the suppression of J c with magnetic field, as considered elsewhere in [14,17,18,21], 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 5 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: 1) as the magnet approaches the centre of the wire from the right-hand side (t = 347 ms in the 2 nd cycle); 2) when the magnet is aligned with the centre of the wire (t = 353 ms); and 3) as the magnet moves away from the wire on its left-hand side (t = 359 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 3, with 1) close to the first negative peak, 2) close to the trough in between the 2 nd and 3 rd (positive) peaks and 3) close to the fourth negative peak. At 1), 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 2), these forward and reverse currents are almost equal and opposite, such that V ≈ 0, and at 3), the reverse situation to 1) 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.

Discussion
As shown in Figures 3-5, 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 5.5) and, where possible, the number of mesh elements in the HTS wire set to 120 x 1 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 10 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,54] to increase the HTS layer thickness from 1 µm to 100 µ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 4) and tighter tolerances settings could not improve this performance dramatically. The spikes seen in the J and E plots in Figure 5 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,55] for the T-A (1D) model used for AC loss calculations. It was recommended in both [45,55] 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,55], this did not impact the calculations of interest significantly (AC loss and voltage, respectively), although the voltage waveforms in Figure 3 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 3D models in particular, the scalar potential formulation introduces 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 selfprogrammed 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,57]. 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 3.1.1) 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 5.4, ran significantly slower when opened and run in version 5.5. 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 5.5, 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. 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 opencircuit 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 equationbased equivalent circuit methods.
Excellent qualitative and quantitative agreement was obtained between all models for the opencircuit 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 10 cycles of PM rotation, was calculated to be -9.41 µV with a standard deviation of 0.34 µ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 10 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 1 µm to 100 µ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.