Overview of H-Formulation: A Versatile Tool for Modeling Electromagnetics in High-Temperature Superconductor Applications

This paper reviews the modeling of high-temperature superconductors (HTS) using the finite-element method (FEM) based on the <inline-formula> <tex-math notation="LaTeX">$\boldsymbol {H}$ </tex-math></inline-formula>-formulation of Maxwell’s equations. This formulation has become the most popular numerical modeling method for simulating the electromagnetic behavior of HTS, especially thanks to the easiness of implementation in the commercial finite-element program COMSOL Multiphysics. Numerous studies prove that the <inline-formula> <tex-math notation="LaTeX">$\boldsymbol {H}$ </tex-math></inline-formula>-formulation is able to simulate a wide scope of HTS topologies, from simple geometries such as HTS tapes and coils, to more complex HTS devices, up to large superconducting magnets. In this paper, we review the basics of the <inline-formula> <tex-math notation="LaTeX">$\boldsymbol {H}$ </tex-math></inline-formula>-formulation, its evolution from 2D to 3D, its application for calculating critical currents and AC losses as well as magnetization of HTS bulks and tape stacks. We also review the use of the <inline-formula> <tex-math notation="LaTeX">$\boldsymbol {H}$ </tex-math></inline-formula>-formulation for large-scale HTS applications, its use to solve multi-physics problems involving electromagnetic-thermal and electromagnetic-mechanical couplings, and its application to study the dynamic resistance of superconductors and flux pumps.


I. INTRODUCTION
In recent years, the electromagnetic and mechanical performance of high-temperature superconductors (HTS) has remarkably increased. The high current density and the possibility of cryogen-free operation make HTS the most suitable candidate for future superconducting applications [1]. Another advantage of HTS is that their price has been gradually decreasing [2]. All these reasons make HTS a technically and economically viable solution for various superconducting applications [3]- [18], and some of such applications have almost reached the market. While the price of HTS is still too high to compete with well-established conventional applications such as power cables and transformers, certain HTS applications are already well established: resistive superconducting fault current limiters are at a pre-commercial stage and HTS magnets have reached record-breaking magnetic fields.
The associate editor coordinating the review of this manuscript and approving it for publication was Derek Abbott .
Meanwhile, the investigation and research of HTS have also developed from simple objects to more complex devices and to large-scale applications. For simulating the behavior of complex HTS systems, analytical solutions are not accurate, and experimental investigations can be challenging, time consuming, and even risky. Therefore, numerical tools are becoming increasingly important to design HTS applications and predict their performance.
For the simulation of the electromagnetic behavior of HTS, several numerical models based on different formulations, such as A−V [19], T− [20], T−A [21], and H-formulation [22]- [24] have been proposed. In 2003, the H-formulation, which takes its name from the single dependent variable of magnetic field intensity H, was first used to calculate macroscopic current and field distributions in superconductors in home-made codes [22], [25] and successively implemented in the commercial finite-element method (FEM) software COMSOL Multiphysics by Hong et al [23] and Brambilla et al [24]. The H-formulation has been widely used for modeling various HTS topologies, with the advantages of accuracy, good convergence and acceptable VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ computing time. Based on the citations of the two papers [23], [24] mentioned above and on presentations at technical conferences, we have identified more than 45 research groups worldwide using the H-formulation for modeling the electromagnetic behavior of superconductors. This paper overviews the previous works on HTS with the H-formulation FEM models, starting from the basics of the H-formulation, and covering its evolution from 2D to 3D, its application for calculating critical currents and AC losses and for investigating the magnetization of HTS bulks and tape stacks. Reviewed here is also the use of the H-formulation for large-scale HTS applications, its use for solving multi-physics problems involving electromagneticthermal and electromagnetic-mechanical couplings, and its application for studying the dynamic resistance of superconductors and flux pumps. A schematic illustration of some of the possible uses of the H-formulation is given in Figure 1. This review can be helpful for researchers to understand the mechanism, development, and scope of applications of the Hformulation FEM model for HTS, which could be beneficial for the future designs of HTS applications.

II. BASICS OF THE H-FORMULATION
From the mathematical point of view, the H-formulation uses the finite-element method to solve Faraday's equation: using the magnetic field H as state variable and a nonlinear resistivity for describing the peculiar electrical characteristics of superconductors: where E 0 is a characteristic electric field (generally 10 −4 V/m), J is the current density, J c is the critical current density and n is a factor indicating the steepness of the transition from the superconducting to the normal state. This powerlaw model for the superconductor material approaches Bean's critical state model as the power law factor n increases toward infinity [26]. The power-law is commonly used because it describes the power dissipation typically occurring at power frequencies. For slower phenomena (e.g. relaxation), the percolation-law could give a more precise description [27], [28]. The non-conducting regions are modelled as materials with high resistivity. A resistivity value of 1 ·m is often used: it is sufficiently large, thus avoiding the flow of current in non-conducting regions, and it does not excessively increase the computation time.
From the magnetic point of view, the superconductor is modelled as a material with relative magnetic permeability µ r = 1. Faraday's equation expressed in terms of magnetic field takes the form: The zero-divergence of the magnetic flux density is guaranteed by the choice of zero initial conditions, as explained in [29]. The current density J is derived from the magnetic field by means the quasi-static approximation of Ampere's law: Equation (3) can be solved by FEM software programs such as COMSOL Multiphysics [24], Matlab [30], FlexPDE [31], home-made FEM codes such as Daryl Maxwell [32], and open-source environments such as GetDP [33].

III. EVOLUTION FROM 2D TO 3D
The early studies based on the H-formulation mostly employed two-dimensional models, simulating the transverse cross section of long superconductors with the current flowing perpendicular to the cross section. In such models, the magnetic field has only two components (in the plane of the analyzed cross section) and the current density has only one component, in the direction perpendicular to the field. Those early models were used to calculate AC losses and were successfully validated against experimental measurements and analytical solutions [23], [24], [34], [35]. The challenge of simulating superconductors with very large width-to-thickness ratio (such as HTS coated conductors) was first tackled by artificially expanding the superconducting layer's thickness [36] and then, for systems characterized by a large number of tapes, with the homogenization method [37]. In addition, Rodriguez-Zermeno et al developed a useful way to reduce the number of degrees of freedom using a structured and elongated rectangular mesh [38], which allowed accurately simulating the superconducting layer with its real thickness and saving computation time.
The H-formulation was then extended to 3D, in order to handle problems of increasing complexity. In general, the equations are the same, and it is assumed that the current density is parallel to the electric field -an assumption that is not necessarily true in three dimensions. In addition, in the majority of cases, it is assumed that the resistivity given by equation (2) is isotropic.
In [39], Zhang et al built 3D H-formulation FEM models for HTS bulks and HTS tapes, and calculated their electromagnetic profiles in order to estimate the AC losses. In [40], a full 3D model of superconducting twisted wires transporting AC electric current was established. In the same work, the 3D model was also used to analyze the influence of material defects in HTS tapes. In [29], Zermeno et al built a full 3D model of a Roebel cable made of 14 strands, and they did a comprehensive investigation on current distributions and AC losses: they found that a significant advantage of the 3D model over its 2D counterpart was the capability of revealing the presence of small highly dissipative regions near the corners of Roebel cable. In [41], a full 3D model of a racetrack HTS coil was built with the homogenization approach and compared to 2D models. Two-dimensional axisymmetric H-formulation FEM models were also developed [42]- [44].

IV. AC LOSS CALCULATION
AC losses can constitute a serious refrigeration burden, limiting the efficiency and increasing the operational cost of superconducting applications. Their estimation is therefore very important. AC losses in HTS tapes and coils can be measured with different experimental techniques, such as the electrical method [45], the magnetic method [46], and the calorimetric (boil off) method [46]. However, the measurements require dedicated experiments, which are often delicate and lengthy. In addition, AC loss measurements become much more challenging with large-scale superconducting applications. Therefore, reliable tools which are able to accurately estimate the AC loss values of superconducting devices are particularly sought for.
In the past decades, several analytical models for the calculation of AC losses in superconductors have been developed. These models provide close-form expressions for the AC losses as a function of few physical and geometrical parameters. Among the most used ones are the expressions by Norris [47] and Brandt [48] for the transport and magnetization losses in individual superconductors with simple geometry (such as infinitely long tapes and wires with elliptical or thin rectangular cross section). A comprehensive review of available analytical models is given in [49]. Analytical models are efficient for computing the losses of simple geometries of HTS, but in general they cannot be used for large-scale superconducting devices. By contrast, numerical models can model geometries and situations of increasing complexity. One of the main advantages of the H-formulation FEM model is its accurate AC loss computation for various HTS topolo-gies under the actions of different current and magnetic field conditions.
Currently, the most widely used HTS tapes are rare-earth barium copper oxide (REBCO) coated conductors. They consist of different layers: the superconducting layer (which is extremely thin, typically 1 µm), the encapsulation, the stabilizer, and the substrate. In general, there are four types of losses that possibly contribute to the total AC loss of a REBCO HTS tape [50]: • hysteresis losses, in the superconducting material; • coupling losses, when tapes are subdivided into filaments; • magnetic losses, if there are magnetic materials; • eddy-current losses, in the metals. In most cases of practical interest, the hysteresis losses in the superconducting material dominate, if the frequency of the AC excitation is not too high (e.g. less than 1 kHz) and there are no magnetic materials or filaments. The H-formulation is primarily used for calculating the hysteresis losses, but the other types of losses can also be calculated [50]. In the H-formulation model, the AC loss can be calculated by integrating the instantaneous power dissipation over the domain of interest (the superconductor or the metal) and taking the average over a cycle [51]- [53]: where T is the time interval of each cycle of AC signal, and is the domain for loss calculation. For the purpose of more accurate loss results, various versions of B-dependent VOLUME 8, 2020

FIGURE 2.
Comparison of the H-formulation FEM model, measurement, and analytical solutions (Norris ellipse and Norris strip), for the analysis of the AC loss from a HTS a single tape and double pancake coil, with transport current at 50 Hz, with data from [56].
critical current models are usually incorporated into the H-formulation [54], [55]. Figure 2 shows results for the AC losses of a single HTS tape and a double pancake HTS coil fabricated with the same material (SCS6050 from SuperPower), with a 50 Hz AC transport current [56]. The AC losses calculated with the H-formulation for both single tape and coil agree with the experiments on a wide range of values (more than four orders of magnitude).
The H-formulation can calculate the losses from the different materials of REBCO HTS tapes, such as the magnetic AC losses in the magnetic substrates, and the eddy-current AC losses in the metals (e.g. copper stabilizer and silver over-layer). Nguyen et al used a modified version of the H-formulation with a non-linear magnetic permeability to calculate the total AC losses in the HTS tapes with magnetic substrates [57]. Shen et al used the H-formulation model to compute the eddy-current AC losses in the copper stabilizer, and the overall AC losses in the HTS tapes (SuperPower, SCS12050 and SF12100). The results were verified by the measurements [58].
The H-formulation can be used to calculate the electromagnetic interaction and the AC losses of tapes when they are assembled in the form of vertical stacks or horizontal arrays. The magnetization AC losses of a vertical stack of HTS tapes were studied in [34], with a good match between the simulation results and the experimental values. Shen et al carried out experiments and simulations for the AC losses in an horizontal array of HTS tapes and predicted the AC losses with increasing numbers of tapes in the array [51].
Several studies have been dedicated to the calculation of AC losses by means of the H-formulation. Hong et al modelled the magnetization losses of a superconducting racetrack coil with the conditions of ripple field together with DC background field or DC transport current [59]. Xu and Grilli investigated the AC losses in a 4-mm HTS tape in the presence of various DC and AC ripple currents [60]. Lahtinen et al simulated ripple field losses under the action of direct current and varying magnetic field, finding a good agreement with experiments in most of the considered cases [61]. In general, care should be taken for evaluating AC ripple losses [62]. The H-formulation uses a power-law resistivity (see equation (2)), which causes a relaxation of the current density distribution associated to the DC component. The value of the cyclic AC ripple losses calculated by means of equation (5) depends on when the calculation is performed along such relaxation. Usually one has to wait until the relaxation of the DC component does not influence the results. This usually requires simulating a large number of cycles.
Nguyen et al carried out an AC loss study on anti-parallel connected HTS tapes, which is a beneficial configuration for superconducting fault current limiter applications [63]. A method to reduce AC losses of HTS tapes is by striating the tape into a large number of filaments and narrow grooves: the efficiency of this method was demonstrated with experiments and confirmed by FEM simulations in [64]. Grilli et al used the H-formulation to carry out the numerical modeling and loss analysis of a twisted stacked tape cable (TSTC) [65]. The possibility of considering the effect of contact resistances on the current distributions in HTS cables was introduced in [66].
The H-formulation FEM model can be a useful tool to simulate various configurations of HTS cables, such as the Roebel cable, conductor-on-round-core (CORC) cable, and other ones. In [67], Grilli and Pardo calculated the AC losses of Roebel cables, and their study provided insights on the losses for the cases of coupled and uncoupled HTS strands. Several AC loss studies on Roebel cables were performed by means of both experiments and simulations based on the H-formulation (both 2D and 3D) with different conditions of transport current and magnetic field [29], [68]- [71]. Majoros et al and Terzioğlu et al investigated the magnetization and transport AC losses of HTS CORC by experiments and simulations [72], [73]. There are other types of HTS cables, such as the twisted stacked tape (TSTC) HTS wires, soldered-stacked-square-twisted (3S-T) HTS cable, quasi-isotropic strands, HTS Cross-Conductor Cables, whose AC loss characteristics have been studied by using the H-formulation [65], [74]- [76].
Quéval et al used the H-formulation to calculate the AC loss of large-scale HTS coated conductor applications. Their work also included the comparison of a full model (simulating all the turns of a large coil) with multi-scale and homogenized models [77]. These last two models -presented in more detail in Section VII -allow dramatically reducing the size of the problem and the computation time, without significant loss of accuracy.
A systematic review of the AC loss computation for HTS using the H-formulation is presented in [78].

V. CRITICAL CURRENT AND DC CHARACTERISTICS
The H-formulation can also be used to estimate the critical current and DC characteristics of HTS tapes and coils. In [79], Zhang et al used the 2D axisymmetric H-formulation to study the magnetic field distribution of an HTS coil, and determined the coil's critical current by considering the weakest turn of the coil. As a result of the analysis, they suggested a coil (end-to-end) voltage of 20 µV/m (instead of the conventional 100 µV/m) as a new criterion for determining a long term and safe operating current. They also found that the presence of an anti-phase background field changes the position of the weakest turn of the coil. In [69], the 2D axisymmetric H-formulation was used to simulate pancake coils made from a Roebel cable and calculate the critical current density of each HTS layer and the magnitude of the magnetic flux density in and around the coil. In [70], Thakur et al studied the voltage-current characteristics and critical currents of superconducting strips and Roebel cables. In [55], Grilli et al proved that dynamic simulations using a current ramp and evaluating the average electric field in the superconductor can be a proper method to determine the critical current: their results were in good agreement with a self-consistent DC model and with measurements.
In [80], Ainslie et al calculated the DC characteristics and critical current of an HTS coil with magnetic flux diverter under different background magnetic fields. In [81], the DC properties were successfully simulated by a 2D axisymmetric H-formulation model, and the anisotropy and non-uniformity of HTS coils were studied. In [82], Hu et al used the measurement of angular dependence of the critical current of HTS tapes, and incorporated the J c (B, θ ) characteristics into the Hformulation FEM model. Their model could accurately match the critical current and angular dependence of a SuperPower HTS tape (SCS4050), and they found that simulations with the two-variable direct interpolation of the J c (B, θ ) data were over 100 times faster than those using a data fitting method. In [83], Liu et al carried out a comprehensive comparison of critical current simulations of a double pancake racetrack coil with over 200 turns: the H-formulation was used as a comparison for other models, which can be computationally more efficient for this kind of calculation: the T-A formulation (which considers the superconductor layer as an infinitely thin object), and the P-model and modified load-line method (which use static simulations). The predictions of the models were similar to each other and in reasonably good agreement with the experiments.

VI. MAGNETIZATION OF BULKS AND STACKS
In the fundamental studies of HTS modeling, the magnetization pattern of a HTS bulk or strip can be well modelled by H-formulation and the results are in good agreement with analytical solutions [23], [24]. In [84], the magnetization curves of an HTS bulk in the presence of different amplitudes and frequencies of magnetic field were simulated in 2D. In [39], Zhang and Coombs used the 3D H-formulation to model a single HTS bulk and a bulk array, and showed that HTS bulk arrays can create a highly uniform magnetic field profile in a reasonably large domain. In [85], Zhang et al examined the pulse magnetization of an HTS bulk and studied the distribution of trapped field using the H-formulation.

In [86], Kapolka et al studied different magnetization behaviors of HTS rectangular bulks and stacks in the presence of inclined fields using different 3D models, including the H-formulation.
In [87], Ainslie et al performed a pulsed field magnetization (PFM) experiment for two c-axis oriented, single grain GdBCO and YBCO bulk superconductors, and used the 3D H-formulation model to investigate the effects of inhomogeneities on the trapped field and maximum temperature increase of the samples. In [88], a novel arrangement of the HTS bulk with magnetizing coils and iron yoke to trap high field was experimentally tested, and a 2D axisymmetric H-formulation model was used to qualitatively reproduce and monitor the magnetization process from the thermal and electromagnetic points of view. In [89], Ainslie and Fujishiro presented an overview of the modeling of HTS bulk magnetization, and proved that both the 3D and 2D axisymmetric H-formulation models are viable tools for various conditions of HTS bulk research.
In [90], Page et al carried out the experiment of pulsed magnetization for HTS stacks using multiple pulses with different temperatures, and used the 2D H-formulation and the real thickness of each layer in the HTS stack to accurately model the process of eddy current and heating effects. In [91], Baskys et al simulated trapped magnetic fields and critical current density of HTS stacked tapes with angular transversal field. In [92], [93], Zou et al studied the influence of the n index and magnetic field constant B 0 value in the Kim's model on the field trapped in HTS bulks by means of pulsed field magnetization. An example of an HTS bulk coupling cooling boundary and thermal insulation in the presence of pulsed magnetic field is presented in Figure 3, which shows the modeling scheme and the domains and the boundary conditions (Figure 3(a)), the temporal evolution of the applied field pulse (Figure 3(b)) and some results at selected instants of the pulse (Figure 3(c)). In the last plot, the influence of the index n on the results for the current density distribution is clearly visible, as are the different results obtained with a purely electromagnetic model (on the left) and a coupled electromagnetic-thermal model (on the right). In [94], the same authors investigated an HTS stack magnetized using pulsed field magnetization by controlling the current of external coils based on H-formulation. Later, Zou et al performed experiments on HTS stack magnetization to validate the electromagnetic-thermal coupled 2D H-formulation model, and the results indicated that using suitable magnetization sequences could efficiently improve the trapped field [95].

VII. LARGE-SCALE HTS APPLICATIONS
The H-formulation can be used to simulate the superconducting parts of large-scale HTS applications, in particular to investigate their electromagnetic performance (e.g. loss analysis) and to study the effect of screening currents.
Shen et al used the H-formulation FEM model to develop superconducting magnets using a Halbach array configuration for transportable medical imaging device called Lorentz Force Electrical Impedance Tomography (LFEIT) [96], and then they used optimization methods to further improve the uniformity of the superconducting magnet [97]. As shown in Figure 4, Shen et al also incorporated the superconducting magnetic properties into the LFEIT system to simulate the magneto-acoustic signal from different biological samples [98]. Xia et al evaluated the electromagnetic performance and calculated the AC losses of a highfield superconducting magnet using the H-formulation with anisotropic bulk homogenization model, with the purpose of getting insights for the superconducting prototype coils of the National High Magnetic Field Laboratory 32 Tesla allsuperconducting magnet [99].
The H-formulation can be a useful tool for modeling superconducting fault current limiters. In [100], a 10 kV resistivetype HTS fault current limiter was investigated by means of experiments and H-formulation calculations. In [101], Jia et al performed a loss and magnetic field analysis of a saturated iron-core superconducting fault current limiter (SISFLC). In addition to the spatial distribution of the AC losses, they studied the relationship between the AC component of the current and the AC losses. In [102], Shen et al investigated the details of the power dissipation of a threephase 35 kV/90 MVA SISFCL. The losses were estimated to be up to the kW level and potentially even higher, depending on the amplitude of the used DC bias current.
Wang et al modelled a hybrid HTS magnet with around 7000 turns using the homogeneous bulk approximation. Such magnet was proposed to be applied in a superconducting magnetic energy storage (SMES) system [103]. Using the same homogeneous bulk approximation, Song et al simulated an HTS three-phase 1 MVA transformer, and compared the loss results with those previously obtained with the minimum magnetic energy variation (MMEV) model, finding a reasonably good agreement [104].
Several studies have been reported on the modeling of HTS motor or generator applications based on the H-formulation. Hu et al used a 3D model to simulate all-superconducting synchronous electric machines, and they performed the analysis on the magnetic field around the HTS coils and proposed a method to reduce the total AC losses [105]. Zhang et al simulated the HTS racetrack armature windings for largescale HTS machines and studied their magnetization and transport current AC losses [106], [107]. For HTS generators, Quéval et al combined the power grid analysis and the loss study of a 10 MW class HTS wind turbine generator [108]. Li et al modelled the armature structures of a fully HTS synchronous generators with loss and electromagnetic analysis [109]. Brambilla et al developed a hybrid A-H formulation to investigate the electromagnetic behavior of HTS electrical machines [110]. The idea behind this approach was the separation of the model of an electrical machine in two parts, where the magnetic field is calculated with the most appropriate formulation: the H-formulation in the part containing the superconductors and the A-formulation in the part containing conventional conductors (and possibly permanent magnets).
The work focused on determining and applying the continuity conditions on the boundary separating the two regions.

VIII. COUPLING WITH THERMAL MODELS
In COMSOL Multiphysics, the H-formulation can be conveniently coupled with other physics modules. A typical thermal model is based on the time-dependent heat transfer equation: where T is the temperature, λ and C v are the thermal conductivity and the volumetric heat capacity. Q ac is the heat dissipation generated by the AC loss which can be calculated using the H-formulation. Roy et al used a thermal coupled H-formulation to model the thermal and electromagnetic behavior of HTS tapes for resistive superconducting fault current limiters [111]. Majoros et al studied the heat generation and magnetization AC losses in a HTS CORC by measurements and simulations based on the 2D thermal-coupled H-formulation [72]. Song et al investigated the AC losses and their impact on the thermal stability of HTS windings in superconducting machines [112]. Huang et al presented a numerical study of the YBCO bulk for maglev application together with loss analysis in [113].
Hu et al built a thermal-coupled 2D model for HTS tapes using homogenization to analyze the effects of the substrate layer on the lightning current performance of HTS tapes. The simulation results were verified by experiments and proved that the maximum temperature always occurred at the top of the homogenized bulk [114].
Zou et al used the electromagnetic-thermal coupled 2D axisymmetric H-formulation bulk models to study the impact of the magnitude and frequency of external AC magnetic fields on the bulk's trapped field profile [115]. As mentioned in section VI, these studies generally incorporated the thermal-coupled model into H-formulation FEM model to achieve more realistic physical characteristics [87]- [89], [94], [95].

IX. COUPLING WITH MECHANICAL MODELS
The mechanical-coupled H-formulation can be used to investigate magnetic forces, levitation and suspension of HTS applications. Ma et al developed an analytical model based on the H-formulation solver coupled with force calculation for the magnetic levitation and suspension of HTS applications, and they validated the analytical solutions with experiments [117], [116]. Liu et al performed measurements of magnetic levitation with an HTS stacks made of 120 REBCO coated conductor tapes, and showed that the 2D H-formulation model was able to accurately predict the levitation performance [118]. Huang et al studied the influence of both intrinsic factors (e.g. critical current density) and external factors (e.g. the moving velocity, and moving path of HTS bulk) on the levitation performance of HTS maglev systems, based on current density distribution and magnetic levitation analysis [119].
Sass et al proved that the H-formulation FEM model can be used to calculate the levitation force on HTS stacks and bulks of superconducting bearing applications and that it is a proper approach to model moving objects [120], and the vertical levitation force in y-direction per unit length can be calculated by: where S T is the surface of all superconducting regions. Grilli et al used H-formulation and the Arbitrary Lagrangian-Eulerian formulations to model the levitation phenomena of superconducting bulks and permanent magnets with a dynamic mesh approach [121], and Zheng et al used a similar method to investigate the electromagnetic characteristics of superconducting maglev devices [122]. Ta and Gao built a 3D electromagnetic-mechanical H-formulation FEM model to study the current density, the magnetic field distribution and the Lorentz force of a Nb 3 Sn superconducting strand [123]. Patel et al used a force-, field-and temperaturedependent H-formulation model to analyze the induced current density distribution of HTS annuli, made from tapes with ferromagnetic substrate. They found that the presence of the ferromagnetic substrate could enhance the stable axial levitation force to values higher than possible with a uniform bulk of infinite J c [124].

X. DYNAMIC RESISTANCE AND FLUX PUMPS
The dynamic resistance and flux pumping characteristics of HTS can be modelled with the H-formulation. Ainslie et al developed a 2D model to compute the dynamic resistance VOLUME 8, 2020 and total loss in an HTS coated conductor transporting a DC electric current and in the presence of background AC magnetic fields. This model helped explain the origin and characteristics of dynamic resistance including the relationship with the applied field and current, the onset of nonlinear flux-flow resistance and the threshold field [125]. As shown in Figure 5, Geng et al used a 2D H-formulation model to simulate a transformer-rectifier flux pump using the mechanism of switching dynamic resistance. Figure 5(a) shows a typical transformer-rectifier circuit: a superconducting charging loop with varying magnetic field and a superconducting load, and Figure 5(b) shows the corresponding schematic for the FEM modeling using H-formulation. Their analysis showed that the model was capable of simulating the flux pumping behaviors and the simulation results had good agreement with the experimental data [126]. Mataira et al used the H-formulation to verify the principle of DC output voltage of an HTS dynamo, and demonstrated the process of large overcritical currents flowing within the HTS stator during certain periods of each dynamo cycle [127].

XI. ADVANTAGES AND DISADVANTAGES
The H-formulation has shown its powerful capability to describe the electromagnetic behavior of a wide range of HTS applications and it is now the standard FEM model used for that purpose. There are several reasons for the popularity of the H-formulation. First and foremost, the H-formulation is easy to implement in the widespread commercial software package COMSOL Multiphysics. This is an important reason, because it opens this kind of simulations to researchers that are not necessarily experts in numerical modeling techniques. Since a few years, the built-in Magnetic Field Formulation (MFH) module in COMSOL Multiphysics makes the use of the H-formulation model even simpler with respect to when equation (3) had to be written in the general Partial Differential Equation (PDE) module. Another advantage of the MFH module with respect to the traditional PDE approach is that it is easier to simulate axisymmetric geometries and consider materials with nonlinear magnetic properties.
In addition, thanks to the widespread use of the Hformulation in COMSOL Multiphysics, it is easier for research groups and the whole community to share FEM models and transfer the knowledge between members of a research group. Some models are publicly available on the website of the HTS Modeling Workgroup [128], which can be useful to learn about problem settings like mesh size, tolerances, etc.
Despite the success, the H-formulation has some important limitations: (i) For large-scale applications involving HTS coated conductors, the T-A formulation [21], [129], which simulates the superconducting layer as an infinitely thin object, is more efficient than the H-formulation: for example, in [21] the T-A formulation was reported to be about 20 times faster than the H-formulation for the simulation of 500 individual HTS tapes.
(ii) The H-formulation simulates the air domains, and the calculation speed is slower than other methods which do not do that [67].
(iii) This puts also some limits on the number of tapes that can be simulated in 2D, which is typically in the range of several hundred [77]. For larger numbers, the homogenization method is a viable solution, as long as the superconductor layer can be considered infinitely thin (no currents induced across the thickness).
(iv) The possibility of using parallel computing with the H-formulation needs to be better addressed, especially for the solution of large 3D problems [130].
As a final remark, one should also note that, from the mathematical point of view, a more appropriate approach to solve this kind of problems would use scalar potentials in the non-conductive regions and vector fields only in the conductive regions. The advantage of these hybrid formulations should become evident, especially for solving 3D problems of increasing complexity. Some initial steps in that direction have recently been taken [131].

XII. CONCLUSION
This article presented a systematic review on the simulation of high temperature superconductors and their applications using the H-formulation FEM model. This model has been and still is primarily used for calculating AC losses, but it can also be used for other purposes, like simulation of DC characteristics and calculation of trapped fields and flux pumps. The simulation of multi-physics problems, particularly electromagnetic-thermal and electromagneticmechanical coupled problems, is facilitated by the implementation in the COMSOL Multiphysics finite-element program. The H-formulation has been also implemented in other programming environments, but they represent a minority.
Numerous studies demonstrated that the H-formulation FEM model is one of the most reliable numerical modeling methods to simulate various HTS applications -from simple topologies like HTS bulks, tapes, coils, to devices and largescale applications -with good agreement with experiments and a reasonable computation speed. This paper reviewed the previous HTS studies based on H-formulation FEM models, starting from the basics of the formulation, and following its evolution towards cases of increasing complexity, both from the geometrical and physical point of view. Although the recently developed T-A formulation is becoming popular for the simulation of HTS coated conductors, the H-formulation is much more flexible and, for the foreseeable future, it is expected to retain a major role for the simulation of superconducting applications.