<figure><img src="image_27.png" /><figcaption>Fig. 6. Transient transmitted electric field.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 7. Shielding effectiveness.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 8. Transient transmitted electric field.</figcaption></figure>

The difference is from the homogenization model to explain the error between the case with interface condition and the case when the sheet is meshed [10], but this difference is acceptable and the model remains valid.

## 6 Thee dimensional cavity

Let's consider a 30 cm × 30 cm × 30 cm cavity with 3 cm × 3 cm aperture located on the front size (Fig. 10). The wall is constituted by two composite sheets. Each sheet is constituted by fibers reinforced by resin. The characteristics of the resin are $ε_1 = 2ε_0$, $μ_0$, $σ_1 = 0.01$ S/m and

<figure><img src="image_30.png" /><figcaption>Fig. 9. Shielding effectiveness.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 10. 3D cavity.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Fig. 11. Homogenization steps.</figcaption></figure>

characteristics of the fibers are $ε_2 = ε_0$, $μ_0$, $σ_2= 10000$ S/m. The thickness of the sheets is 6 mm. The fibers are oriented in the same direction. The sheets are assembled such that the fibers belonging to two adjacent sheets are oriented perpendicularly (Fig. 11) [13]. Each heterogeneous sheet is homogenized separately (Figs. 11.1.a and 11.1.b). The obtained homogenized sheets are associated (Fig. 11.2). A homogenization procedure is performed again to obtain an equivalent homogeneous wall (Fig. 11.3).

The homogeneous wall is anisotropic and his properties depends on the polarization of the electromagnetic field. They are defined by the orientation of the normal vector (Wiener model [5] or dynamic homogenization model).
<figure><img src="image_33.png" /><figcaption>Fig. 12. Electric field in the center of the cavity.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 13. Electric field in the center of the cavity.</figcaption></figure>

In this study, the normal vectors of the walls are oriented along the canonical basis ($e_x, e_y, e_z$):

- where the normal vectors are oriented along $e_x$ and $e_z$, the effective parameter are obtained with the Wiener model, the effective permittivity is $\epsilon = 2.38\epsilon_0$, and the effective conductivity is $\sigma = 990$ S/m;

- where the normal vectors are oriented along $e_y$, the effective parameter are presented in Figures 4 and 5.

The cavity is illuminated by a Gaussian pulse, $E = (0, E_y, 0)$ and $H = (H_x, 0, 0)$, which excites frequencies up to 1.5 GHz. The magnitude of the Gaussian pulse is 1. The simulations are performed on a computer with 32 CPU (2 GHz Pentium Xeon, Intel, Santa Clara, California, USA). The size of the mesh is about 15 000 elements. An explicit Runge-Kutta 44 scheme is used. The CFL time step is $4 \times 10^{-14}$ s. The time cost of one is simulation is about 5 days.

The electric field in the center of the cavity is plotted in Figure 12, and compared to the case where the walls of the cavity are perfect electric conductor (PEC) and the case walls of the cavity are conductive $\sigma = 5$ S/m (the thickness of the walls is the same than composite walls). The electric field remains important with composite walls because the conductivity of the walls oriented along $e_y$ is low. The attenuation of the electric field in the case of

conductive walls is more important than of the case with PEC walls. A FFT is applied and the shielding effectiveness is evaluated and presented in Figure 13. A resonance frequency is excited in the PEC case and is attenuated in the cases of conductive and composite walls. The shielding effectiveness is less important than the PEC case due to the transmission of the walls.

This simulation demonstrates the efficiency of the approach for studying a cavity with heterogeneous walls. If the walls of the cavity were meshed, the computational cost would be strongly increased due to the presence of microstructures, and the simulation would become cumbersome.

## 7 Conclusion

An interface condition to avoid the mesh of heterogeneous sheets in 3D geometries is proposed in a DG method. This interface condition is issued from analytical solution of 1D Maxwell equations in frequency domain. It allows a reduction of the computational time cost. This interface condition was validated on homogenization cases. A study of a 3D cavity was presented to illustrate the efficiency of the approach.

## References

1. K.M. Mitzner, IEEE Trans. Antennas Propag. **16**, 706 (1968)

2. F. Bocquet, L. Pichon, A. Razek, IEEE Trans. Magn. **34**, 2791 (1998)

3. M.S. Sarto, IEEE Trans. Electromagn. Compat. **41**, 298 (1999)

4. M. Boubekeur, A. Kameni, L. Bernard, A. Modave, L. Pichon, IEEE Trans. Magn. **50**, 493 (2014)

5. O. Wiener, Abhandlungen der Kniglich Schsischen Gesellschaft der Wissenschaften **32**, 604 (1912)

6. A. Sihvola, IEEE Electromagnetic Waves Series 47, (1999)

7. G.W. Milton, The Theory of Composites (Cambridge University Press, 2002)

8. C.L. Holloway, M.S. Sarto, M. Johansson, IEEE Trans. Electromagn. Compat. **47**, 833 (2005)

9. V. Preault, R. Corcolle, L. Daniel, L. Pichon, IEEE Trans. Electromagn. Compat. **55**, 1178 (2013)

10. V. Preault, R. Corcolle, L. Daniel, L. Pichon, J. Appl. Phys. **115**, 154904 (2014)

11. A. Kameni Ntichi, A. Modave, M. Boubekeur, V. Préault, L. Pichon, C. Geuzaine, Eur. Phys. J. Appl. Phys. **64**, 24508 (2013)

12. T. Mori, R. Tanaka, Acta Metall. Mater. **21**, 597 (1973)

13. G. Wasselynck, D. Trichet, J. Fouladgar, IEEE Trans. Magn. **49**, 1825 (2013)
## 2.7 Diagnostic de réseau filaire

Les réseaux filaires occupent une place importante dans les systèmes tels que la voiture, l'avion ou encore le train. La fiabilité de ces réseaux repose sur la qualité des outils et techniques de maintenance employés, notamment ceux de la maintenance prédictive car un dysfonctionnement dans le réseau peut conduire à des accidents mettant en péril la sécurité des biens et des personnes. La détection de petites dégradations de câbles avant que celles-ci ne deviennent plus importantes est d'un intérêt majeur pour la fiabilité de ces réseaux. Cette problématique autour des défauts naissants est notamment étudiée au laboratoire GeePs à travers des travaux de thèse de Mostapha SMAIL [104] et a bénéficié de plusieurs financements ANR dont SODDA. Ma contribution sur ce thème a consister en la mise en place et à l'utilisation des modèles numériques 3D pour la caractérisation de ces défauts naissants [67, RI.4, RI.7].

La technique la plus couramment utilisée pour le diagnostic des réseaux filaires est la réflectométrie. Celle-ci consiste à injecter un signal dans le réseau et à étudier les signaux réfléchis (Fig.2.4) dont la présence dénote de l'existence d'un ou plusieurs défauts [101, 103]. Les caractéristiques telles que la nature et la position des défauts détectés sont obtenus après l'analyse de la forme et du déphasage des signaux réfléchis [42, 91].

Dans les travaux de modélisation, les reflectogrammes utilisés pour caractériser les défauts

<figure><img src="image_35.png" /><figcaption>FIGURE 2.4 – Schéma de principe de la réflectométrie.</figcaption></figure>

d'un réseau filaire sont très généralement issus de la résolution des équations des télégraphistes. C'est un système d'équations mono-dimensionnelles dont les inconnues sont la tension V(z, t) et le courant I(z, t). Il est issu de la représentation quadripolaire électrique d'une ligne de transmission en fonction des éléments suivants : R (Résistance), L (Inductance), C (Capacité), G (Conductance), et est établi en utilisant la loi des noeuds et la loi des mailles sur un tronçon de câble de longueur $\Delta z$ (Fig.) [90].

$$ egin{cases} \frac{\partial V}{\partial z} + L \frac{\partial I}{\partial t} + RI(z, t) = 0 \\ \frac{\partial I}{\partial z} + L \frac{\partial V}{\partial t} + GV(z, t) = 0 \end{cases} \quad (2.34) $$
Pour déterminer les paramètres électriques RLCG nécessaires au système des télégraphistes,

<figure><img src="image_36.png" /><figcaption>FIGURE 2.5 – Représentation RLCG d'une ligne de transmission.</figcaption></figure>

on utilise généralement deux types d'approches qui sont définies uniquement sur la section de
la ligne étudiée :

- Une approche analytique : elle donne les paramètres électriques en fonction des caractéristiques géométriques du câble [106]. Dans le cas de câbles co-axiaux sans pertes, l'inductance et la capacité sont données par : $L = (\mu/2\pi) \ln (r_a/r_b)$, $C = 2\pi\epsilon/\ln (r_a/r_b)$, où $r_a$ est le rayon du conducteur intérieur et $r_b$ le rayon du conducteur extérieur.

- Une approche numérique 2D : elle consiste à résoudre équation de poisson $\nabla \cdot (\epsilon \nabla V) = 0$ [65]. Le potentiel V solution de ce problème statique est utilisé pour le calcul des paramètres électriques de la ligne.

Malheureusement, les paramètres électriques calculés par ces deux approches ne tiennent pas compte de l'aspect tri-dimensionnelle du câble et des défauts. L'un des objectifs de nos travaux a été de mettre en évidence l'apport de la modélisation 3D sur la détermination de ces paramètres électriques RLCG [RI.7].

### 2.7.1 Modèle de réflectométrie pour la caractérisation des défauts non-francs

Pour mener nos travaux sur l'étude des défauts non-francs, un modèle numérique tri-dimensionnel basé sur la résolution des équations de Maxwell en domaine temporel a été mis en place et ses résultats ont été confrontés aux mesures expérimentales.

Le modèle numérique proposé est défini sur un tronçon de câble contenant le défaut à étudier. Un champ électrique incident est appliqué à l'une des extrémités du câble et les signaux réfléchis et transmis sont enregistrés sur des plans situés en avant et en arrière du défaut. A partir de ces signaux temporels, les coefficients de réflexion et de transmission du défaut sont calculés. Le modèle expérimental quant à lui a consisté à fabriquer des défauts le plus proche possible de ceux modélisés, et à mesurer leurs coefficients de réflexion grâce à un analyseur de réseau.
Contrairement aux réflectogrammes numériques, les réflectogrammes expérimentaux sont bruités par des signaux parasites non désirés. Des techniques permettant d'isoler la réponse du défaut sont mise en oeuvre et rendent les résultats expérimentaux comparables à ceux de la simulation [46].

Dans les développements réalisés, on s'est particulièrement intéressé aux trois type de défauts non-francs suivant :

- Le défaut d'écrasement qui traduit une variation de la section du câble sur une longueur $L_f$ (Fig.2.6).

- Le défaut de blindage qui traduit l'arrachement du blindage sur un angle $\theta_f$ et une longueur $L_f$ (Fig.2.7).

- Le défaut d'épaulement qui traduit l'arrachement du blindage et du di-électrique sur une profondeur $p_f$ et une longueur $L_f$ (Fig.2.8).

<figure><img src="image_37.png" /><figcaption>FIGURE 2.6 – Défaut d'écrasement</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>FIGURE 2.7 – Défaut de blindage</figcaption></figure>

Le cas du défaut d'épaulement est présenté dans l'article ci-après **Time domain reflectometry model : analysis and characterization of a chafing defect in a coaxial cable** publié dans le journal *The European Physical Journal, Applied Physics*. La comparaison aux résultats expérimentaux réalisée montre la validité de l'approche numérique mise en place. Ces résultats sont exploités pour l'évaluation des paramètres électriques équivalents dans le but de proposer une représentation quadripolaire plus fine prenant en compte les aspects 3D.
# Time domain reflectometry model: analysis and characterization of a chafing defect in a coaxial cable*

Abelin Kameni*, Florent Loete, and Lionel Pichon

GeePs – Group of electrical engineering – Paris, UMR CNRS 8507, CentraleSupélec, Université Paris-Sud, Université Paris-Saclay, Sorbonne Université, 11 rue Joliot-Curie, Plateau de Moulon, 91192 Gif-sur-Yvette Cedex, France

Received: 22 December 2017 / Received in final form: 7 June 2018 / Accepted: 11 June 2018

**Abstract.** This paper presents experimental and numerical studies of a chafing soft defect realized by partially milling coaxial cables. The approach is based on the time domain reflectometry technique. The numerical model consists in solving Maxwell's equations while an incident Gaussian pulse is injected on the faulty line. The experimental time domain measurements are performed with a vector network analyzer. To get the experimental results comparable to the numerical ones, a process to denoise the measured impulse responses is proposed. The reflection coefficients obtained are compared to those given by a classical approach based on a chain matrix model to show the impact of 3D numerical modeling in studying soft faults.

## 1 Introduction

One of the crucial issues in domains such as aerospace and automotive is to keep electrical networks undamaged by detecting and repairing the defects that can appear on the wires [1]. The hard faults which are characterized as open-circuit or short-circuit have been well studied and are efficiently characterized by traditional reflectometry techniques [2]. The soft faults which are characterized by partial damage to the wire have been far less studied because the generated reflections are very small and hard to detect [3]. To improve the analysis and diagnosis of wired networks, it is crucial to have numerical models whose results are comparable to those obtained from measurements.

Modeling a wiring system generally consists in computing the voltage potential on the faulty section by solving a 2D static Poisson equation and using the Gauss theorem to determine the distributed electrical parameters: the resistance *R*, the inductance *L*, the capacitance *C*, and the conductance *G* (RLCG). These parameters are used in a longitudinal model such as telegrapher's equations or chain matrix model to compute the reflection coefficient [4]. Unfortunately, this modeling approach does not take into account the three-dimensional aspect of wave propagation in the wire. In our previous works, we have proposed how to compute the equivalent electrical parameters of shielding soft defects through 3D time domain numerical modeling and use frequency domain reflectometry measurements for validation [5]. Unfortunately, in these works, because of the considered fault lengths and the cable type, the description of reflection

coefficient is limited on the growing part of the first lobe. This does not allow to clearly determine the frequency signature of the shielding defect.

In this paper, we use time domain reflectometry measurement and a 3D time domain numerical simulation to describe evolution of the frequency signature of soft fault defects. The damage is a partial chafing of the shield and the insulator on a coaxial cable RG58. The fault lengths and the frequency band are chosen sufficiently large to have several lobes on the reflection coefficient. In this work, we also show how to treat experimental time domain measurements to get results which can be compared to those obtained from numerical modeling. Because of very weak contribution of the soft defect on the response of the line, the proposed process combines different existing techniques such as de-embedding [6] and clean algorithm [7], in order to remove the unwanted signals on the measured impulse responses. The experimental and numerical reflection coefficients due to the soft defect are compared to that obtained from a chain matrix model to show benefit of 3D modeling.

## 2 Numerical modeling

A 3D numerical modeling consists in solving the time domain Maxwell's equations to compute the reflection coefficient of the faulty section. A nodal discontinuous Galerkin method is adopted for the spatial discretization [8]. Its discontinuous aspect allows to easily discretize objects of different sizes or shapes, and provides a better representation of discontinuous properties. These kind of methods are well adapted for parallel computing and high order spatial elements are easily used to improve the accuracy of the results.

<footnote>* Contribution to the topical issue “Numelec 2017”, edited by Adel Razek.</footnote>

<footnote>* e-mail: abelin.kameni@lgep.supelec.fr</footnote>
## 2.1 Numerical method

Let **E**, **H**, and **J** represent, respectively, the electric field, the magnetic field, and the current density. The time domain Maxwell's equations form a system (1) of 6 unknowns that are components of **E** and **H**:

$$ \left\{ \begin{array}{l} \varepsilon \partial_t \mathbf{E} - \nabla \times \mathbf{H} = -\mathbf{J} \\ \mu \partial_t \mathbf{H} + \nabla \times \mathbf{E} = 0 \end{array} \right. , \qquad (1) $$

where $\varepsilon$ is the permittivity of the medium, $\mu$ its permeability. In a conductive medium, $\mathbf{J} = \sigma \mathbf{E}$, with $\sigma$ the conductivity.

The discontinuous Galerkin method is introduced for solving the conservative form of partial differential equations. This method consists in discretizing the variational formulation of (1) on each mesh element $T$ of the domain $\Omega = \cup T$:

$$ \left\{ \begin{array}{l} \int_T \varepsilon \partial_t \mathbf{E} \varphi - \int_T \mathbf{H} \times \nabla \varphi - \int_{\partial T} (\mathbf{n} \times \mathbf{H}) \varphi = - \int_T \sigma \mathbf{E} \varphi \\ \int_T \mu \partial_t \mathbf{H} \psi + \int_T \mathbf{E} \times \nabla \psi + \int_{\partial T} (\mathbf{n} \times \mathbf{E}) \psi = 0 \end{array} \right. , \qquad (2) $$

where $\varphi$ and $\psi$ are test functions. In each $T$, a finite element method is applied and a mapping technique is employed to facilitate the use of high order mesh elements. The interface terms $(\mathbf{n} \times \mathbf{H})$ and $(\mathbf{n} \times \mathbf{E})$ are replaced by numerical flux expressions $(\mathbf{n} \times \mathbf{H})^{\text{num}}$ and $(\mathbf{n} \times \mathbf{E})^{\text{num}}$ as in a finite volume method. Different formulations of the flux expressions exist [9]. These following expressions (3) resulting in different numerical schemes are implemented. For $\alpha = 0$, centered fluxes are obtained and numerical schemes are dispersive. For $\alpha = 1$, upwind fluxes are obtained and numerical schemes are dissipative. We have

$$ \left\{ \begin{array}{l} (\mathbf{n} \times \mathbf{H})^{\text{num}} = \mathbf{n} \times \frac{\{\mathbf{ZH}\}}{\{\mathbf{Z}\}} - \alpha \left( \mathbf{n} \times \frac{\{\mathbf{n} \times \{\mathbf{E}\}\}}{\{\mathbf{Z}\}} \right) \\ (\mathbf{n} \times \mathbf{E})^{\text{num}} = \mathbf{n} \times \frac{\{\mathbf{YE}\}}{\{\mathbf{Y}\}} + \alpha \left( \mathbf{n} \times \frac{\{\mathbf{n} \times \{\mathbf{H}\}\}}{\{\mathbf{Y}\}} \right) \end{array} \right. , \qquad (3) $$

where $Z = \sqrt{\mu/\varepsilon}$, $Y = 1/Z$, $\|\mathbf{u}\| = \frac{u^+ - u^-}{2}$, and $\{u\} = \frac{u^+ + u^-}{2}$. The superscript “-” denotes the values for fields in the current element, while “+” is for the adjacent element. In this work, the upwind fluxes are considered, the time integration is ensured by an explicit four stages Runge-Kutta method (RK44) and the spatial discretization is performed with the third order tetrahedral mesh elements to improve accuracy of the results.

## 2.2 Numerical results

The coaxial cable model is designed by four concentric cylinders as shown in Figure 1a. At the center, the copper core diameter is 0.43 mm. The first layer of diameter 2.9 mm is the inner dielectric insulator whose relative permittivity is $\varepsilon_{r_i} = 2.3$. The second layer of diameter 3.2 mm is the woven copper shield whose conductivity is $\sigma = 6.10^7 \text{ S/m}^{-1}$. The third layer of diameter 3.6 mm is an

<figure><img src="image_39.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Fig. 2. Impulse responses of a chafing soft fault of length $L_f = \{5, 7.5\}$ cm.</figcaption></figure>

outer plastic jacket whose relative permittivity is $\varepsilon_{r_j} = 2.1$. The soft fault presented in Figure 1b is a partial chafing of the jacket, the shield, and the insulator on a length $L_f$, and the damaged depth of the dielectric insulator is noted as $p_f = 1$ mm. Since a part of the field is radiated outside, a surrounding cylindrical vacuum volume is added with perfectly matched layers to operate like free space as in experimental measurements [10].

The simulation consists in injecting a Gaussian pulse of frequency band [0, 7] GHz in the cable. The incident radial electric field is given by: $E_r^{\text{inc}}(t) = e^{-a(t-t_0)^2}$ with $a = 3.10^{20}$ and $t_0 = 0.3$ ns. The reflected field is recorded on an observation plane $S$ located in front of the soft fault. It is computed as

$$ E_r^{\text{refl}} = S^{-1} \int E_r dS. \qquad (4) $$

Two lengths $L_f = \{5, 7.5\}$ cm are considered and the impulse responses versus the position $z = 2vt$, with $v = c/\sqrt{2.3}$, are presented in Figure 2. For each case,
<figure><img src="image_42.png" /><figcaption>**Fig. 3.** A two port a vector network analyzer (VNA) is connected to a faulty line through the fixtures A and B.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Fig. 4.** The faulty line includes a chafing soft defect realized by partially removing by milling the shielding and the isolator.</figcaption></figure>

the first reflection is positive and the second is negative. This is due to the higher wave velocity value in the defective zone. The distance between the positive and negative peaks is equal to the length of the damaged part. The amplitudes of the negative pulses are lower than those of the positive pulses. This is due to the radiated field outside the coaxial cable. When $L_f$ increases the amplitude of the second reflection decreases because the field is radiated outside along the defective zone.

# 3 Experimental approach

In this section we will focus on the comparison between the numerical simulation and the experimental results on chafing defects for diagnosis purposes. The early diagnosis of such defects is very important since they correspond to electrical wire damaged by abrasion against sharp edges caused by poor design, bad or broken holding of the cable. This results in a chafing of the insulator as modeled in the previous section and ultimately, in an extreme case where the insulator is completely removed over an angular portion, can lead to severe system failures.

## 3.1 Experimental setup

The experimental setup is shown in Figure 3. A faulty line is connected to a vector network analyzer through two 50Ω coaxial cables equipped with N connectors. The A and B fixtures are N to BNC coaxial adaptors. Since the SOLT (Short Open Load Thru) calibration is realized at the extremities of those cables, the position of the reference plane will be at the end of the cable connected to port 1. The scattering parameters are measured with an Agilent N9914 Fieldfox over a [0, 4]GHz range.

<figure><img src="image_44.png" /><figcaption>**Fig. 5.** Original and de-embedded impulse responses of a defective line under test. The line is a 60 cm 50Ω RG58 coaxial cable including a 7.5 cm chafing defect.</figcaption></figure>

In this work, the faulty line including a chafing defect (Fig. 4) is realized by partially milling the shielding and the insulator on a limited portion of a 60 cm, 50Ω RG58 coaxial cable, as shown in Figures 1 and 3. Two lengths $L_f$ = {5, 7.5} cm with a depth $d$ = 1 mm were realized.

## 3.2 Experimental results

A common problem, when trying to characterize a soft defect on an electrical line by reflectometry, is that the fixtures at each extremity of the line also introduce insertion or return losses. In that case, the measured impulse responses also contain the signatures of the fixtures. Although very weak, the fixtures contribution cannot be neglected when dealing with soft defects since the perturbation they create on the line are of the same order of magnitude. In this work the two-window time gating method presented in reference [6] is used to remove the contribution of the fixtures. It is based on the measurement of the four scattering parameters of the network and a selective gating of the mismatches caused by each fixture. Consequently the full de-embedding of the fixtures effect and the quantitative characterization of the signature of defects can be achieved. The result of the de-embedding process on the defective line is presented in Figure 5 and one can see that the contribution of the fixtures are mostly removed. Despite this technique, it remains a non-negligible noise level that leads on ripplings of the reflection coefficient in frequency domain.

To avoid these ripples, a clean algorithm procedure [7] is applied in order to remove the noise and the unwanted inhomogeneities on the time domain impulse response. In this work, the shape of the impulse response at one interface of the defect can not be considered as a *sinc* function. As a matter of fact, the leakage of the field with frequency will affect the reflection coefficient and consequently the shape of the impulse response. The interface impulse response used in the clean algorithm procedure was taken as the main lobe of the first reflection (Fig. 6).

Figure 7 shows the comparison of the reflection coefficients between the experimental measurements and the numerical model for $L_f$ = {5, 7.5} cm. The error over
<figure><img src="image_45.png" /><figcaption>**Fig. 6.** Result of the de-embedding and clean algorithm procedure on the impulse responses of the defective line under test with a chafing defect of length $L_f = 7.5$ cm.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Fig. 7.** Comparison of reflection coefficient between the experimental measurements, the numerical DG model and a simple chain matrix model for $L_f = 7.5$ cm and $d = 1$ mm.</figcaption></figure>

frequency is around 5% below 1.5 GHz and increases with frequency since many experimental parameters can have an influence on the losses. The reflection coefficient obtained from a chain matrix model, that consists in a transmission line model with the three sections cascaded, is also added for $L_f = 7.5$ cm. This chain matrix model does not take into account the losses due to the leakage in the outer medium, and can account for the model of a crushing defect where the shielding would not have been damaged. The periodicity of the lobes with frequency is in very good agreement but this result more importantly clearly emphasizes the effect of the losses on the frequency response. Consequently, the specific signature of the chafing defect evidenced here is very important since it

allows us to discriminate a simple crushing defect (no leakage, no short circuit risk) from a chafing defect potentially more severe. The length of the defect can be estimated through the measurement of lobes periodicity and the identification of the nature of the defect through the recognition of the attenuation pattern by cross correlation analysis.

## 4 Conclusion

A numerical model, as well as an experimental validation, of chafing soft defects on an electrical line was presented. The model is based on the implementation of the discontinous Galerkin method and the defects created by milling coaxial cables were characterized by time and frequency domain reflectometry. A de-embedding and clean algorithm procedure was also used in order to denoise, extract, and quantitatively identify the signature of the soft defect. The method is very efficient and, for the first time to our knowledge, a frequency signature allowing to discriminate a crushing from a chafing defect was evidenced. This result can be a powerful tool for the preventive diagnosis and maintenance of potentially severe soft defects.

## References

1. F. Auzanneau, Prog. Electromagn. Res. B **49**, 253 (2013)

2. M. Franchet, N. Ravot, O. Picon, IEEE Electromagn. Compat. Mag. **2**, 54 (2013)

3. F. Loete, Q. Zhang, M. Sorine, IEEE Trans. Antennas Propag. **63**, 2532 (2015)

4. E.J. Lundquist, J.R. Nagel, S. Wu, B. Jones, C. Furse, IEEE Sens. J. **13**, 1172 (2013)

5. A. Manet, A. Kameni, F. Loete, G. Genoulaz, L. Pichon, O. Picon, IEEE Trans. Electromagn. Compat. **59**, 533 (2017)

6. G. Gronau, I. Wolff, IEEE Trans. Microw. Theory Tech. **37**, 479 (1989)

7. C. Buccella, M. Feliziani, G. Manzi, IEEE Trans. Electromagn. Compat. **46**, 597 (2007)

8. A. Kameni, A. Modave, M. Boubekeur, V. Preault, L. Pichon, C. Geuzaine, Eur. Phys. J. Appl. Phys. **64**, 24508 (2013)

9. J.S. Hesthaven, T. Warburton, J. Comput. Phys. **181**, 186 (2002)

10. A. Modave, A. Kameni, J. Lambretchs, E. Delhez, L. Pichon, C. Geuzaine, Eur. Phys. J. Appl. Phys. **64**, 24502 (2013)

**Open Access** This article is distributed under the terms of the Creative Commons Attribution License https://creativecommons.org/licenses/by/4.0 which permits unrestricted use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

**Cite this article as:** Abelin Kameni, Florent Loete, Lionel Pichon, Time domain reflectometry model: analysis and characterization of a chafing defect in a coaxial cable, Eur. Phys. J. Appl. Phys. **83**, 30601 (2018)
<figure><img src="image_47.png" /><figcaption>FIGURE 2.8 – Défaut de blindage et de di-électrique</figcaption></figure>

### 2.7.2 Evaluation des paramètres électriques équivalent

Pour mettre en évidence l'apport de la modélisation 3D, on a comparé le coefficient de réflexion issu du modèle numérique à celui obtenu d'un modèle de ligne de transmission utilisant les paramètres (RLCG)an analytiques. Dans le cas d'un défaut d'écrasement sur une ligne sans pertes, les travaux de thèse de Hossein Manesh (Anthony Manet) ont montrés qu'il y a un décalage en fréquence sur les coefficients de réflexion Fig.2.9. Ce décalage est dû à une différence sur la valeur de l'élément capacitif.

L'observation d'une cartographie du champ électrique du modèle numérique 3D, montre la

<figure><img src="image_48.png" /><figcaption>FIGURE 2.9 – Comparaison des coefficients de réflexion : (g) modèle numérique 3D et modèle des télégraphistes avec les paramètres (RLCG)an, (d) modèle numérique 3D et modèle des télégraphistes avec les paramètres (RLCG)num</figcaption></figure>

présence d'effets capacitifs dus à la propagation d'une composante $E_z$ au niveau des arêtes de changements de section du câble (Fig.2.10). A partir de l'évaluation des paramètres électriques
<figure><img src="image_49.png" /><figcaption>FIGURE 2.10 – Cartographie du champ électrique</figcaption></figure>

$(RLCG)_{num}$, une solution d'amélioration de la représentation électrique quadripolaire a pu être proposé une amélioration. Celle-ci consiste en l'ajout d'éléments capacitifs en entrée et en sortie du quadripôle représentant la zone défectueuse (Fig.2.11). Cette solution proposé a permis d'éliminer le décalage en fréquence observé précédemment (Fig.2.9).

<figure><img src="image_50.png" /><figcaption>FIGURE 2.11 – (g) Représentation quadripolaire classique, (d) Représentation quadripolaire avec ajout d'éléments capacitifs</figcaption></figure>

Lorsque le défaut comporte une perte de blindage (défaut de blindage et épaulement), le modèle 3D montre qu'il y a un rayonnement vers l'extérieur du câble. Ce rayonnement conduit à l'atténuation des ondes à la traversée de la zone défectueuse. Lorsque la fréquence augmente, celui-ci devient plus important ce qui se traduit par la diminution de l'amplitude des lobes du coefficient de réflexion. Ces travaux de thèse on permis de mettre en évidence que en plus des
effets capacitifs dus aux changements de section du câble, il y a des effets dissipatifs. Pour prendre en compte ces effets dissipatifs, la solution proposé dans ces travaux est d'ajouter dans la représentation électrique de la zone défectueuse, un élément résistif qui est mis en série avec l'inductance. Dans le tableau présenté ci-après, on récapitule les modèles électriques équivalent améliorés des trois types de défauts étudiés. L'article ci-après **Equivalent Circuit Model of Soft Shield Defects in Coaxial Cables Using Numerical Modeling** publié dans le journal *IEEE Transactions On Electromagnetic Compatibility* présente la méthodologie et les résultats obtenus pour le cas du défaut de blindage avec détérioration de la couche diélectrique.
<table><thead><tr><th rowspan="3">Type du câble<br>Type du défaut</th><th colspan="4">Câble coaxial RG-58 (50 Ω, 66% cc)</th><th colspan="4">Câble coaxial W/D (50 Ω, 83% cc)</th><th colspan="4">Câble coaxial W/N (50 Ω, 81% cc)</th></tr><tr><th colspan="4">Ecrasement</th><th colspan="4">Blindage</th><th colspan="4">Isolant</th></tr><tr><th colspan="4">Pf en % de Df (Lf = 2.5 cm)</th><th colspan="4">Bf en degré (Lf = 5 mm)</th><th colspan="4">Pf en % de Df (Lf = 5 mm)</th></tr></thead><tbody><tr><td rowspan="2">Sévérité du défaut</td><td colspan="2">Sain</td><td>10%</td><td>14%</td><td>38%</td><td colspan="2">Sain</td><td>50°</td><td>90°</td><td>130°</td><td colspan="2">Sain</td><td>40%</td><td>90%</td></tr><tr><td colspan="2">L' (uH/m)</td><td>0.248</td><td>0.234</td><td>0.228</td><td>0.185</td><td colspan="2">Ls (uH/m)</td><td>0.189</td><td>0.193</td><td>0.204</td><td colspan="2">Ls (uH/m)</td><td>0.193</td><td>0.210</td><td>0.228</td></tr><tr><td rowspan="2">Paramètres électriques</td><td colspan="2">C' (pF/m)</td><td>100</td><td>106</td><td>109</td><td>134</td><td colspan="2">Cs (pF/m)</td><td>81.52</td><td>80.02</td><td>77.74</td><td colspan="2">Cs (pF/m)</td><td>84.55</td><td>79.92</td><td>73.38</td></tr><tr><td colspan="2">C<sub>d</sub> (pF)</td><td>0.0</td><td>0.0010</td><td>0.0016</td><td>0.009</td><td colspan="2">C<sub>d</sub> (pF)</td><td>0.0</td><td>0.007</td><td>0.015</td><td colspan="2">C<sub>d</sub> (pF)</td><td>0.0</td><td>0.042</td><td>0.108</td></tr><tr><td rowspan="2">f<sub>1</sub> = 3 GHz<br>f<sub>2</sub> = 6 GHz</td><td colspan="2">R<sub>s</sub> (Ω/m)</td><td colspan="3">0.0</td><td colspan="2">f<sub>1</sub><sup>2</sup></td><td>0.0</td><td>0.037</td><td>1.361</td><td colspan="2">f<sub>1</sub><sup>2</sup></td><td>0.0</td><td>3.546</td><td>28.91</td></tr><tr><td colspan="2">f<sub>2</sub><sup>3</sup></td><td colspan="3">0.0</td><td colspan="2">f<sub>2</sub><sup>3</sup></td><td>0.0</td><td>0.353</td><td>10.44</td><td colspan="2">f<sub>2</sub><sup>3</sup></td><td>0.0</td><td>11.31</td><td>93.99</td></tr><tr><td>Module du coefficient<br>de la réflexion</td><td colspan="2">|Γmax|</td><td colspan="3">0.0</td><td colspan="2">|Γmax|</td><td>0.0</td><td>0.021</td><td>0.071</td><td colspan="2">|Γmax|</td><td>0.0</td><td>0.061</td><td>0.141</td></tr><tr><td>Impédance du défaut</td><td colspan="2">%Zd<sup>4</sup></td><td colspan="3">0.0%</td><td colspan="2">%Zd<sup>4</sup></td><td>0.0%</td><td>+2.25%</td><td>+2.76%</td><td colspan="2">%Zd<sup>4</sup></td><td>0.0%</td><td>+2.63%</td><td>+3.67%</td></tr><tr><td>forme du signal<br>réfléchi</td><td colspan="4"></td><td colspan="4"></td><td colspan="4"></td><td colspan="4"></td></tr><tr><td>Schéma équivalent</td><td colspan="4"></td><td colspan="4"></td><td colspan="4"></td><td colspan="4"></td></tr><tr><td>Modèle de la<br>géométrie du défaut</td><td colspan="4"></td><td colspan="4"></td><td colspan="4"></td><td colspan="4"></td></tr></tbody></table>

1. Lf : est la longueur du défaut

2. f1 : fréquence du travail = 3 GHz

3. f2 : fréquence du travail = 6 GHz

4. Variation nominal en %

<caption>Tableau 5-2: Etude paramétrique de l'évolution des caractéristique du défaut par rapport au type et leur degré de la sévérité.</caption>
# Equivalent Circuit Model of Soft Shield Defects in Coaxial Cables Using Numerical Modeling

Anthony Manet, Abelin Kameni, Florent Loete, Jérôme Genoulaz, Lionel Pichon, and Odile Picon

**Abstract**—This paper presents a three-dimensional modeling study to investigate the effects of one type of soft fault problem in electrical fault diagnosis, the impacts of partial degradation in the coaxial cable shielding or other shielded lines. At first, a time domain nodal discontinuous Galerkin method is used to solve Maxwell's equations and simulate reflectometry responses of soft faults. Then, the numerical simulation results for soft defects in cable shields are compared with experimental measurements for validation purposes. Finally, the shield defects that can arise on coaxial cables are modeled by an equivalent circuit and compared to illustrate their impacts on the electrical characteristic parameters as a function of their severity.

**Index Terms**—Cable shielding, coaxial cables, discontinuous Galerkin, equivalent circuits, fault diagnosis, reflectometry, simulation, soft fault.

## I. INTRODUCTION

THE fault-detection feature is certainly a very important aspect of wire health monitoring and an important process required in electrical wiring system operation. It has a great influence on security and quality. In transmission line networks, this feature is needed for the early identification of the faulty line in order to anticipate severe consequences that can be costly and even catastrophic.

Two main categories of defects can occur on transmission lines [1]: hard faults like short or open circuits resulting from a significant change in line impedance and the soft faults that include partial degradation such as insulation damage, natural aging, water infiltration, etc, leading to very small impedance changes along transmission lines. These defects can be studied by the reflectometry technique, which is based on electrical signal injection into the network from one or both extremities of a line followed by an analysis of the reflected signals resulting from impedance discontinuities [2], [3]. Whereas hard faults are easy to find, soft faults are difficult to detect because they do not significantly affect the received signal in

traditional reflectometry measurements. These soft faults have been less studied and they are more difficult to identify because of their low-level reflection signatures until they become hard faults [4]–[6]. In order to anticipate the appearance of severe faults that are initially caused by soft fault degradation, a detection feature is needed to provide timely identification of the faulted line.

The objective of this paper is to carry out a soft fault model in a coaxial line system to analyze its electrical effects caused by a discontinuity in the shield. In a previous study [7], the authors have presented a preliminary discussion on a numerical study of a soft defect by computing their reflectometry responses. This paper aims at extending the applicability of such a technique to provide an equivalent circuit for the shield defect and show the influence of fault severity on electrical parameters. The numerical results presented in this paper are obtained with the discontinuous Galerkin method. This method has been recently introduced in the modeling of electromagnetic compatibility problems and has proved to provide a powerful tool for wave problems [8]. The unstructured mesh allows a flexible discretization of complex geometries and provides adequate consideration of discontinuous physical properties. This method is well adapted for parallel computing because the generated matrices are block diagonals. The easy use of high order mesh elements allows for increasing the accuracy of the results.

## II. NUMERICAL MODELING OF A SOFT DEFECT

The aim of this section is the numerical simulation of electromagnetic fields in the case of a soft defect, in order to describe the behavior of the reflection line response in the presence of soft faults. We refer to a lossless coaxial cable transmission line of characteristic impedance $Z_0$, with the damaged area between two interfaces separated with a discontinuity in the shield of length $L_f$.

The application example was modeled in three dimensions by means of the discontinuous Galerkin method applied to the soft fault problem to numerically solve the time-dependent Maxwell's equations [8], [9], and to compute different parameters such as the time-domain reflection response, the reflection coefficient, or the impedance of the faulty section.

The discontinuous Galerkin methods are introduced for solving the conservative form of partial differential equations. This method consists of discretizing the variational formulation on each element of the computational domain. In each element, a finite-element method is applied, and the mapping technique is used to facilitate the use of high-order mesh elements. The interface terms are replaced by numerical flux expressions as in a finite-volume method [10], [11].

In this study, the coaxial cable model provides dimensionally accurate and detailed 3-D models of shielded defects (see Fig. 1).

Manuscript received March 11, 2016; accepted April 14, 2016. Date of publication October 13, 2016; date of current version December 6, 2016. This work was supported by the Association Nationale de la Recherche under Project "Soft Defects Diagnosis in Wired Networks."

A. Manet and O. Picon are with the Laboratory of Electronics, Systems of Communications and Microsystems (EA2552), Université Paris-Est Marne-la-Vallee, Marne-la-Vallée 77454, France (e-mail: antony.manet@gmail.com; odile.picon@univ-mlv.fr).

A. Kameni, F. Loete, and L. Pichon are with the Group of Electrical Engineering of Paris, CentraleSupélec, Université Paris-Sud, UPMC, Gif-sur-Yvette 91192, France (e-mail: abelin.kameni@supelec.fr; florent.loete@lgep.supelec.fr; lionel.pichon@geeps.centralesupelec.fr).

J. Genoulaz is with the Safran Electrical and Power, Department of Research and Technology, Electrical Wiring Interconnection System Divison (Eurasia), Blagnac 31700, France (e-mail: jerome.genoulaz@safrangroup.com).

Color versions of one or more of the figures in this letter are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TEMC.2016.2612719
<figure><img src="image_51.png" /><figcaption>Fig. 1. Coaxial cable model: cross-sectional model (a) with different layers and (b) with the soft fault that is a shield and jacket cutaway.</figcaption></figure>

<table><caption>TABLE I
COAXIAL CABLE MODEL DEFINITION</caption><thead><tr><th>Parameter</th><th>Material</th><th>Diameter (mm)</th><th>Electrical properties</th></tr></thead><tbody><tr><td>Core</td><td>Solid silver-plated copper</td><td>2.33</td><td>PEC</td></tr><tr><td>Insulation</td><td>PTFE</td><td>6.0</td><td>v<sub>p</sub> = 0.83 nominal</td></tr><tr><td>Shield</td><td>Silver-plated copper</td><td>7.10</td><td>PEC</td></tr><tr><td>Jacket</td><td>FEP</td><td>7.70</td><td>ε<sub>r</sub> = 2.1</td></tr></tbody></table>

to describe the effects of the soft fault in a realistic coaxial line system.

Three concentric cylinders define the geometry: 1) outer plastic jacket; 2) double-braided shield; and 3) inner dielectric insulator made from aerated fluoropolymer insulation (PTFE). Fig. 1 shows the mesh used to simulate the coaxial cable with dimensions listed in Table I. At the center, the solid silver plated copper core is considered a perfect electrical conductor (PEC) boundary condition. Note that the double-braided shield is considered as a PEC due to its excellent shielding coverage and meshed because it is located between two dielectric materials.

The considered cable length is $L = 5$ cm and a soft fault occurs at $x_f = 2.25$ cm from the beginning of the line, with lengths of 5, 10, and, 20 mm defects on the shield, at various widths (angle cutaways “θ”) of 50°, 90°, and 130°, as shown in Fig. 1. The example given in this paper pertains to vary low loss coaxial cable with signal loss 105 dB/ 100 m up to 8 GHz. Due to the very short length of our cable model, the losses are neglected in the present analysis. Since a part of the field radiates outside the cable due to the discontinuous shield, a surrounding cylindrical air volume truncated with absorbing Silver–Muller boundary conditions is added to the computational domain and operates like free space as in experimental measurements.

The simulation consists in injecting a Gaussian pulse $U_i$ of the unit amplitude of frequency band [0; 6 GHz] in the cable. The reflected and transmitted fields are recorded on the observation planes located in front of and behind the soft fault. Fig. 2 illustrates the time-domain response of the faulty line as a function of the fault angle values “θ” and the length of the defect “$L_f$”.

For both cases in Fig. 2, the shape of the reflected signal exhibits a positive and negative reflection due to the sign of the reflection coefficient at the discontinuity interfaces at $x = x_1$ and $x = x_2$ (see Fig. 1). The shape and amplitude of the time domain response due to the defect depend on the nature and severity of the fault. Moreover, as shown in these figures, the amplitude of the reflected waves is generally very small. The time domain result is significantly lower than the incident signal

<figure><img src="image_52.png" /><figcaption>Fig. 2. Signal reflected for different simulated cases versus: (a) fault angle with $L_f = 10$ mm and (b) fault length with $\theta = 90^\circ$.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>Fig. 3. Reflection coefficients results versus: (a) fault angle with $L_f = 10$ mm and (b) fault length with $\theta = 90^\circ$.</figcaption></figure>

(here about 5% of $U_i$), and therefore, the reflected pulse cannot be detected easily by means of the measurement devices. As shown in Fig. 2(a), the reflected signal can change significantly due to the radiation leakage that increases with the angle of the shielding defect “θ”. In the same way, Fig. 2(b) shows that the amplitude of the reflected signal is increased with the length of the defect. The observed phase shift between the reflected waves of the simulated cases shown in Fig. 2(b) is due to the length of fault. For the fault of 20 mm length, the second reflection appears later than for the defect with smaller length.

The time-domain responses in Fig. 2 allow us to compute the reflection coefficients $\Gamma(f)$ in the frequency domain [12]. Fig. 3 shows the reflection coefficient $\Gamma(f)$ in a frequency range from 0 to 6 GHz for the two aforementioned cases above, against the fault angle and length, respectively.
<figure><img src="image_54.png" /><figcaption>Fig. 4. Shield damage at 8.4 m on a 25 m coaxial cable, with different lengths.</figcaption></figure>

### III. EXPERIMENTAL VALIDATION

In this section, the experimental measurements that were performed with discontinuities in the coaxial cable shield are presented. The measurement was carried out using an aircraft coaxial cable EN 04604 WD ($Z_0 = 50 \pm 2 \Omega$ at 200 MHz). This type of coaxial cable is designed for high-frequency radio communication applications in aeronautic environments with frequency operation up to 8 GHz. As shown in Fig. 4, the faults were made of 130° shield cutaway sections, with lengths of 1 and 2 cm. The line terminals were attached with 50-Ω N series RF coaxial cable connectors. These N connectors used to join coaxial cables were designed to carry signals at radio frequencies up to 18 GHz. The reflected signals from the cable under test were measured by means of a VNA (Rohde & Schwarz ZNB8). In the case of a measurement with $\Delta f = 100$ kHz and a span of 6 GHz, the time domain range is sufficient to prevent the responses from overlapping or aliasing for reflection measurements for a 25 m cable length. The reflection and transmission coefficients (S-parameters) for a healthy cable and a shield damage located at 8.4 m on a 25 m cable were measured by means of a full 2-port VNA, with improved TRL calibration to acquire highly accurate measurements. The frequency-domain response is a mix of all the signals reflected by the discontinuities present in the cable such as mismatch effects and thus it is difficult to estimate the contribution of the sole fault. However, using both damaged and undamaged results, a time domain filter can be applied to remove unwanted signal components, in order to compute the reflection coefficient corresponding to the reflections initiated by the fault. Then, the measured data filtered and gated in the time domain are transformed back to the frequency domain to obtain the reflection coefficient of the soft fault ($\Gamma_{Gate}$) [13].

As discussed above, $\Gamma_{Gate}$ is measured at the beginning of the line placed at 8.4 m from the fault position and, consequently, the measured reflection response will be lower than the simulation response of a small section of the cable, because of the loss of the cable and connector. In order to improve the accuracy of the reflection values between the simulation and experimental measurements, we need to compensate the loss for the measured reflection coefficient. In the case of the defect placed halfway along the test cable, the attenuation due to the cable and connector losses for reflection at the halfway point is very nearly equal to the $S_{21}$ transmission loss in the undamaged case. Then, the reflection-gated response loss can be compensated by dividing the gated response by original transmission coefficient $S_{21}$ of

<figure><img src="image_55.png" /><figcaption>Fig. 5. Measurement and simulation reflection coefficient results for shield damage on the coaxial cable (fault angle $\theta = 90^{\circ}$).</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>Fig. 6. Equivalent circuit model at the faulty section.</figcaption></figure>

the undamaged case [13]:

$$ \Gamma_{GateComp}(f) = \frac{\Gamma_{Gate}(f)}{S_{21}(f)}. \qquad (1) $$

In our case, where the defect is not placed halfway along the test cable, the reflection gated and compensated response can be computed by taking the average between two values obtained from (1) at both line extremities. An example of the measured reflection coefficient gated and compensated ($\Gamma_{GateComp}$) is shown in Fig. 5 for the test case with the lengths of 1 and 2 cm defect on shield. In this figure, the results are compared with the numerical simulation results obtained by means of the discontinuous Galerkin method and good agreement is achieved.

### IV. EQUIVALENT CIRCUIT

Circuit models are useful to predict system responses to different types of fault behaviors. This permits the representation of the discontinuity by simple electrical element combinations and, consequently, analysis of its behavior in terms of classical RLCG electrical parameters. As a matter of fact, depending on the nature of the defect, one can expect one or several parameters to be altered.

The interruption of the shield provides an impedance discontinuity $Z_d$, which can be due to radiation and induction effect changes. There are also shunt admittance elements $Y_d$, usually in the form of capacitance elements, which also arise from the change in outer conductor radius [14]. Fig. 6 shows the equivalent circuit model of the cable system at the faulty section, assuming that the electric charge density on the faulty region is uniform and neglecting the fringing fields of the edges of the shield discontinuity.
<figure><img src="image_57.png" /><figcaption>Fig. 7. Comparison of characteristic per-unit-length distributed parameters for the coaxial cable in the undamaged and damaged conditions with $L_f = 5$ mm: (a) resistance, (b) inductance, and (c) capacitance.</figcaption></figure>

The values of the circuit elements containing two shunt elements and one series element can be determined from the chain matrix [14] calculated from the S-matrix. The S-parameters used in this section were calculated from the simulation results presented in Section II.

The chain parameter matrix of the faulty region between two interfaces at $x = x_1$ and $x = x_2$ (see Fig. 1) can be calculated from the overall S-matrix calculated at the beginning and end of the line, and then shifted from the measurement reference planes to the faulty region location. In representing the soft defect in this way, the effect of the fault in the coaxial cable shield can be modeled by the distributed parameters. The application example is related to the coaxial cable with the inner dielectric insulator made from an aerated PTFE whose dielectric loss is lower than PTFE with 1.5E-4 at 10 GHz, and that the conductance G can be neglected.

Fig. 7 shows the comparison of characteristic per-unit-length distributed parameters, *R*, *L*, *C*, for the undamaged lossless coaxial cable and damaged with 5 mm long shield cutaways of 50° and 90°. Fig. 7(a) presents the per-unit-length resistance *R* simulated for the undamaged and damaged model. The loss due to the radiation of electromagnetic waves from discontinuity increases with an increase in frequency.

From the coaxial cable characteristics definition listed in Table I, the per-unit-length inductance and capacitance can be calculated directly from analytical expressions: $L' = 0.189 \mu H/m$; $C' = 81.5 pF/m$. In Fig. 7(b) and (c), the simulated results are compared with the calculated per-unit-length inductance and capacitance in the undamaged and damaged model separately. As expected, the shielding defect modifies the capacitance and inductance of the line, because the presence of the discontinuity in the shield affects the propagation characteristics inside the cable. Fig. 7(b) and (c) represents small

changes to line inductance and capacitance in the case of damage. These may be due to small negative lumped capacitance and positive lumped inductance changes on the line at the place of the defect.

## V. CONCLUSION

In this paper, we presented a soft shield defect model on a coaxial cable to analyze and estimate the reflection response.

The numerical simulation results, performed with a time domain nodal discontinuous Galerkin method, are validated by comparison with experimental measurements, and have shown good agreement over a frequency range that widely covers the band usually used for the diagnosis of faults in electric transmission lines.

The shield defect was modeled from reflectometry responses, in terms of the distributed transmission line parameters. A soft fault could result in variation of distributed parameters along the transmission line. The sensitivity analysis about electrical circuit parameter is derived to infer the nature of this kind of defect in terms of changes to per-unit-length resistance, inductance, and capacitance at the place of the defect, as a function of frequency, and could provide more useful information for diagnosis.

## REFERENCES

[1] L. Griffiths, R. Parakh, C. Furse, and B. Baker, "The invisible fray: A critical analysis of the use of reflectometry for fray location," *IEEE Sensors J.*, vol. 6, no. 3, pp. 697-706, Jun. 2006.

[2] P. Tsai, C. Lo, Y. C. Chung, and C. Furse, "Mixed-signal reflectometer for location of faults on aging wiring," *IEEE Sensors J.*, vol. 5, no. 6, pp. 1479-1482, Dec. 2005.

[3] C. Furse et al., "Frequency-domain reflectometry for on-board testing of aging aircraft wiring," *IEEE Trans. Electromagn. Compat.*, vol. 45, no. 2, pp. 306-315, May 2003.

[4] G. Cerri et al., "Fault locatio of shielded cables: Electromagnetic modelling and improved measurement data processing," *Proc. IEE-Sci., Meas. Technol.*, vol. 152, no. 5, pp. 217-226, Sep. 9, 2005.

[5] L. Abboud, A. Cozza, and L. Pichon, "A matched-pulse approach for soft-fault detection in complex wire networks," *IEEE Trans. Instrum. Meas.*, vol. 61, no. 6, pp. 1719-1732, Jun. 2012.

[6] M. Franchet, N. Ravot, and O. Picon, "Soft fault detection in cables using the cluster time-frequency domain reflectometry," *IEEE Electromagn. Compat. Mag.*, vol. 2, no. 1, pp. 54-69, Jan.-Mar. 2013.

[7] H. Manesh et al., "Experimental analysis and modelling of coaxial transmission lines with soft shield defects," in *Proc. 2015 IEEE Int. Symp. Electromagn. Compat.*, Dresden, 2015, pp. 1553-1558.

[8] L. Pebernet et al., "Discontinuous Galerkin method applied to electromagnetic compatibility problems: Introduction of thin wire and thin resistive material models," *IET Sci., Meas. Technol.*, vol. 2, no. 6, pp. 395-401, Nov. 2008.

[9] J. Hesthaven and T. Warburton, "Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell's equations," *J. Comput. Phys.*, vol. 181, no. 1, pp. 186-221, 2002.

[10] A. Kameni, B. Seny, and L. Pichon, "Multirate technique for explicit discontinuous Galerkin computations of time-domain Maxwell equations on complex geometries," *IEEE Trans. Magn.*, vol. 52, no. 3, pp. 1-4, Mar. 2016.

[11] A. Kameni et al., "Evaluation of shielding effectiveness of composite wall with a time domain discontinuous Galerkin method," *Eur. Phys. J. Appl. Phys.*, vol. 64, no. 2, pp. 24508-24511, 2013.

[12] D. M. Pozar, *Microwave Engineering*. New York, NY, USA: Wesley, 2005.

[13] J. P. Dunsmore, *"Handbook of Microwave Component Measurements: With Advanced VNA Techniques."* Hoboken, NJ, USA: Wiley, 2012.

[14] F. M. Tesche, M. Ianoz, and T. Karlsson, *EMC Analysis Methods and Computational Models*. New York, NY, USA: Wiley, 1997.
# Chapitre 3

## Calcul des champs induits dans les matériaux supraconducteurs

### Sommaire

<table><tr><td>3.1</td><td>Introduction</td><td>86</td></tr><tr><td>3.2</td><td>Modélisation numérique des supraconducteurs</td><td>88</td></tr><tr><td>3.2.1</td><td>Loi de comportement électrique des supraconducteurs HTC</td><td>88</td></tr><tr><td>3.2.2</td><td>Formulations pour le calcul de champs électromagnétiques dans les supraconducteurs HTC</td><td>88</td></tr><tr><td>3.2.3</td><td>Linéarisation de la loi constitutive</td><td>90</td></tr><tr><td>3.3</td><td>Méthodes avec flux pour le calcul de champs dans les supraconducteurs</td><td>91</td></tr><tr><td>3.3.1</td><td>Méthode des volumes finis</td><td>92</td></tr><tr><td>3.3.2</td><td>Méthode de Galerkin discontinue</td><td>98</td></tr><tr><td>3.4</td><td>Modélisation des câbles supraconducteurs torsadés</td><td>114</td></tr><tr><td>3.4.1</td><td>Modèle équivalent des filaments torsadés</td><td>115</td></tr><tr><td>3.4.2</td><td>Réduction du modèle équivalent</td><td>120</td></tr></table>
### 3.1 Introduction

Avec l'évolution démographique, les besoins technologiques permettant de satisfaire les populations ont explosé. Dans les pays industrialisés on a besoin de fabriquer un nombre sans cesse croissant de dispositifs et systèmes pour les domaines tel que, le transport, les communications ou encore la santé. Malheureusement, on s'est progressivement rendu compte de la diminution des stocks de matières premières notamment celle extraites de la planète. Face à la pénurie programmée de certaines d'entre elles à cause de la raréfaction, et des difficultés d'extraction, on se tourne progressivement vers d'autres matières premières dont l'utilisation dans les systèmes passe par l'innovation technologique. Parmi ces matières premières on peut citer, les sources d'énergie telles que le vent, l'eau ou encore le soleil, et les nouveaux matériaux tels que les composites, les piézo-électriques ou encore les supraconducteurs. Dans le cas des supraconducteurs, ils ont vocation à remplacer les matériaux conducteurs, comme par exemple le cuivre, sur les aspects de transport de courant et production de champ magnétique. En effet ces matériaux présentent comme avantages, de conduire du courant avec très peu de pertes Joule, et d'être fortement diamagnétique [113].

Les premières observations de la supraconductivité ont été faites au début du 20ième siècle par Kamerlingh ONNES. L'état supraconducteur s'est manifesté dans certains métaux par la disparition de la résistivité en dessous d'une certaine température appelée température critique. Durant les décennies suivant cette découverte, les températures critiques qui étaient de quelques dizaines de Kelvin n'ont cessé d'augmenter pour dépasser la centaine de Kelvins [14]. C'est cette augmentation de la température critique qui a conduit à l'essor des travaux de recherches sur ce domaine, car les perspectives d'utilisation sont plus importantes. L'innovation technologique atour de ces matériaux s'opère depuis plusieurs années, et on assiste au développement de dispositifs utilisant cette classe de supraconducteurs dits à haute température critique (HTC). On peut citer en exemple, les câbles de transport d'énergie, les limiteurs de courant, les moteurs supraconducteurs ou encore les appareils IRM avec des aimants supraconducteurs [58, 100]. Ces réalisations ont été rendues possibles grâce aux travaux de caractérisation qui ont été mis en place pour décrire les interactions de ces matériaux avec les champs électromagnétiques [44, 62]. Ceux-ci ont notamment permis de comprendre la phénoménologie de l'état supraconducteur et de proposer des lois constitutives décrivant son comportement thermoelectromagnétique [5, 15, 24].

L'introduction des matériaux supraconducteurs dans les dispositifs du génie électrique passe par la mise en place de modèles numériques permettant de décrire l'évolution des champs électromagnétiques dans l'environnement concerné, afin d'apporter des solutions pour le dimensionnement des dispositifs. Ces travaux de modélisation sont d'autant plus importants car l'expérimentation autour de ces matériaux est contrainte par les dispositifs cryogéniques nécessaires au
refroidissement du matériau en dessous de la température critique où se manifeste la supraconductivité. C'est sur cet aspect que se situent mes travaux sur les matériaux supraconducteurs. Ils portent notamment sur le développement d'outils numériques pour décrire l'interaction de ces matériaux avec les champs magnétiques. Cette interaction se manifeste par des courants induits dans le matériau lorsque celui-ci est soumis à une sollicitation extérieure comme par exemple un champ magnétique ou un courant. Les outils numériques que j'ai proposés sont principalement dédiés au calcul des champs dans ces matériaux dans le but d'évaluer les pertes et l'aimantation [56]. Ce calcul de champs se fait via la résolution des équations de Maxwell couplées aux lois constitutives du matériau.

Avec l'intérêt croissant autour de ces matériaux durant les dernières décennies, de nombreux modèles analytiques et numériques ont été proposés pour évaluer les grandeurs induites dans ces matériaux. Les premiers modèles ont été dédiés à l'étude des supraconducteurs à basse température critique dont la loi de comportement est donnée par le modèle de Bean, et ils ont consisté en la résolution de systèmes différentiels quasi-linéaires [23, 17, 72]. Avec la découverte des supraconducteurs à haute température critique qui sont caractérisés par une loi puissance, les modèles mis en place nécessitent la résolution de problèmes différentiels fortement non-linéaires. Des approches analytiques ont été proposées pour des configurations simples 1D et 2D [RI.25, 70, 73], mais plus généralement on a recours à des méthodes numériques. Les méthodes intégrales ont été implémentées pour résoudre ces équations pour des configurations 2D [20, 102], mais à cause des matrices pleines, elles ont été abandonnées au profit de la méthode des éléments finis.

C'est autour de la méthode des éléments finis que se concentrent les travaux de modélisation des supraconducteurs à haute température critique [19, 45, 78, 126]. De nombreux modèles 2D et 3D ont été proposés et des logiciels commerciaux tels que COMSOL Multiphysics et FLUX3D disposent de modules spécifiques pour la modélisation des matériaux supraconducteurs. Bien que cette approche soit efficace pour la plupart des configurations, il en existe néanmoins pour lesquels, la convergence n'est pas garantie, ou qui nécessitent des développements supplémentaires pour maintenir son efficacité. C'est le cas lorsque l'exposant de la loi puissance est élevé, il faut faire face à des oscillations numériques, ce qui nécessite de réduire considérablement le pas de temps utilisé et d'ajouter un terme résiduel constant dans la loi puissance [45]. Lorsque la configuration étudiée conduit à un très grand nombre d'inconnues, cette méthode est moins efficace à cause des opérations d'inversions matricielles très couteuses en espace mémoire. Des techniques de compression matricielle ou la mise en place du calcul parallèle sont des solutions envisagées, mais qui demandent d'importants développements supplémentaires [6, 13, 16]. C'est pour ces raisons que la mise en place d'autres approches numériques offre des solutions alternatives à la méthode des éléments finis. C'est ce que j'ai notamment réalisé en proposant, des modèles de calculs de champs et de pertes dans les supraconducteurs utilisant, la
méthode des volumes finis, la méthode de Galerkin discontinue, et des techniques de réductions de modèle pour traiter les configurations conduisant à de grands systèmes comme les câbles supraconducteurs. Dans les sections suivantes, je vous présente les principales formulations utilisées pour la modélisation des supraconducteurs à haute température critique et mes principales contributions sur ce sujet.

## 3.2 Modélisation numérique des supraconducteurs

### 3.2.1 Loi de comportement électrique des supraconducteurs HTC

Une grande partie des travaux de modélisation des supraconducteurs à haute température critique s'appuie sur une loi constitutive décrivant le comportement électrique de ces matériaux. Celle-ci introduit, une résistivité non linéaire qui dépend de la densité de courant induite telle que $E = \rho(\mathbf{J})\mathbf{J}$ (3.1a), ou une conductivité non linéaire qui dépend du champ électrique tel que $\mathbf{J} = \sigma(\mathbf{E})\mathbf{E}$ (3.1b) :

$$ \rho(\mathbf{J}) = \frac{E_c}{J_c(\mathbf{B}, \theta)} \left\| \frac{\mathbf{J}}{J_c(\mathbf{B}, \theta)} \right\|^{n(\mathbf{B},\theta)-1} \quad (3.1a) $$

$$ \sigma(\mathbf{E}) = \frac{J_c(\mathbf{B}, \theta)}{E_c} \left\| \frac{\mathbf{E}}{E_c} \right\|^{\frac{1}{n(\mathbf{B}, \theta)}-1} \quad (3.1b) $$

avec, $\mathbf{B}$ l'induction magnétique, $\theta$ la température, $J_c(\mathbf{B}, \theta)$ la densité de courant critique, $E_c$ une grandeur constante communément fixée à $E_c = 1\,\mu\text{V/cm}$ et $n(\mathbf{B}, \theta)$ un paramètre caractérisant les phénomènes d'activation thermique dans le matériau [15]. Lorsque ce dernier vaut $n = 1$, on est dans le cas d'un matériau à l'état normal résistif.

### 3.2.2 Formulations pour le calcul de champs électromagnétiques dans les supraconducteurs HTC

Lorsque les conditions qu'imposent les grandeurs critiques sont satisfaites, si on soumet un matériau supraconducteur à un champ magnétique ou lui injecte un courant, il se développe en son sein des courants induits. La description de ces courants induits et leur évaluation est nécessaire le dimensionnement des dispositifs utilisant ces matériaux. L'évaluation numérique de ces courants se fait via la résolution des équations de Maxwell auxquelles on associe la loi de comportement électrique du matériau. Pour la loi magnétique qui relie l'induction magnétique $\mathbf{B}$ au champ magnétique $\mathbf{H}$, on utilise généralement celle du vide $\mathbf{B} = \mu_0\mathbf{H}$, avec $\mu_0$ la perméabilité du vide.

Différentes formulations de ce problème peuvent être mise en place. Les principales formulations
tions rencontrées dans la littérature sont les suivantes :

1. Formulation en **H** :

C'est la formulation la plus communément utilisée pour la modélisation 3D des matériaux supraconducteurs [19, 45, 126]. Le principal avantage vient du fait que la condition de champ magnétique imposé est directement traitée par une condition de Dirichlet non homogène. Cette formulation nécessite l'ajout d'un volume d'air autour du matériau afin de mieux prendre en compte les effets démagnétisants. Elle s'obtient en utilisant la forme résistive de la loi constitutive $\mathbf{E} = \rho(\mathbf{J})\mathbf{J}$ qui devient $\mathbf{E} = \rho(\mathbf{J})\nabla \times \mathbf{H}$. Elle s'écrit :

$$ \mu_0 \frac{\partial \mathbf{H}}{\partial t} - \nabla \times (\rho(\mathbf{J})\nabla \times \mathbf{H}) = \mathbf{0} \qquad (3.2) $$

2. Formulation en **E** :

Cette formulation est très peu utilisée pour la modélisation des matériaux supraconducteurs. Dans nos travaux nous avons montré qu'elle convient pour traiter les configurations sans effets démagnétisants comme par exemple un supraconducteur transportant un courant. Dans ces cas, l'ajout d'un volume d'air n'est pas nécessaire et la condition de courant injecté est transformée en une condition de Neumann non homogène qui est imposée sur directement les bords du matériau [RI.27]. Pour les configurations 3D avec effets démagnétisants, nous avons montré que cette formulation peut être utilisée dans le cas de géométries simples dont on connaît les coefficients démagnétisants [RI.17]. Cette formulation s'établit en utilisant la forme conductive de la loi constitutive $\mathbf{J} = \sigma(\mathbf{E})\mathbf{E}$. Elle s'obtient en appliquant l'opérateur rotationnel sur l'équation de Maxwell-Faraday :

$$ \mu_0 \frac{\partial[\sigma(\mathbf{E})\mathbf{E}]}{\partial t} - \nabla \times \nabla \times \mathbf{E} = \mathbf{0} \qquad (3.3) $$

3. Formulation en **J** :

C'est la formulation qui devrait être résolue lorsqu'on s'intéresse à la densité de courant induite dans ce matériau. Malheureusement, celle-ci n'est presque jamais utilisée car la mise en place des conditions aux bords est un obstacle considérable. Cette formulation s'obtient en introduisant la forme résistive de la loi constitutive dans la formulation en $\mathbf{E}$ (3.3) :

$$ \mu_0 \frac{\partial \mathbf{J}}{\partial t} - \nabla \times \nabla \times (\rho(\mathbf{J})\mathbf{J}) = \mathbf{0} \qquad (3.4) $$

4. Formulation en **A** – **V** :

C'est l'une des formulations les plus utilisées pour la modélisation des supraconducteurs [1, 89, 105, 116]. Elle présente l'avantage de pouvoir être facilement couplée aux équations de circuits pour le calcul des pertes. Elle nécessite comme la formulation en H
l'ajout d'un volume d'air pour la prise en compte des effets démagnétisants. Le principal inconvénient est qu'il y a une inconnue supplémentaire par rapport à la formulation en H. Pour construire cette formulation on utilise les relations exprimant les champs électronique **E** et magnétique **B** en fonction du potentiel scalaire **V** et du potentiel vecteur **A** :
$$ \mathbf{B} = \nabla \times \mathbf{A}, \mathbf{E} = -\frac{\partial \mathbf{A}}{\partial t} + \nabla \mathbf{V}. $$
Elle s'obtient, avec la forme conductive de la loi constitutive et l'hypothèse selon laquelle les courants ne sortent pas du matériau $\nabla \cdot \mathbf{J} = 0$. Cette formulation est donnée par le système suivant :

$$ \left\{ \begin{array}{l} \sigma(\mathbf{E}) \left( -\frac{\partial \mathbf{A}}{\partial t} + \nabla \mathbf{V} \right) - \nabla \times \nabla \times \mathbf{A} = 0 \\ \\ \nabla \cdot \left[ \sigma(\mathbf{E}) \left( -\frac{\partial \mathbf{A}}{\partial t} + \nabla \mathbf{V} \right) \right] = 0 \end{array} \right. \qquad (3.5) $$

5. Formulation en $\mathbf{T} - \phi$ :

Cette formulation est très proche de la formulation en **H**, et est utilisée dans quelques travaux de modélisation [107, 117]. Elle s'établit avec la définition des grandeurs **T**, $\phi$, **J** = $\nabla \times \mathbf{T}$ et $\mathbf{H} = \mathbf{T} + \nabla\phi$. Elle s'obtient avec la forme resistive de la loi constitutive et l'hypothèse d'absence de charges magnétique $\nabla \cdot \mathbf{H} = 0$.

$$ \left\{ \begin{array}{l} \mu_0 \frac{\partial[\mathbf{T} - \nabla\phi]}{\partial t} - \nabla \times (\rho(\mathbf{J})\nabla \times \mathbf{T}) = 0 \\ \\ \nabla \cdot \mathbf{T} - \Delta\phi = 0 \end{array} \right. \qquad (3.6) $$

Chacune de ces formulations donne lieu à un problème différentiel très non-linéaire. Cette non-linéarité est portée, par l'opérateur spatial lorsque la forme résistive $\mathbf{E} = \rho\mathbf{J}$ de la loi constitutive est utilisée (formulations **H** et $\mathbf{T} - \phi$), et par l'opérateur temporel lorsque c'est la forme conductive $\mathbf{J} = \sigma\mathbf{E}$ de la loi constitutive qui est utilisée (formulations $\mathbf{E}$ et $\mathbf{A} - \mathbf{V}$). Les modèles numériques développés pour le calcul des champs induits dans les supraconducteurs doivent faire face aux nombreux problèmes de convergence introduits par ces non linéarités. Dans mes développements numériques réalisés autour des formulations $\mathbf{E}$ et $\mathbf{H}$, je propose d'utiliser une linéarisation de la loi constitutive pour améliorer la convergence.

### 3.2.3 Linéarisation de la loi constitutive

Dans nos développements, nous avons travaillé avec une loi constitutive dont les principaux paramètres sont constants, c'est à dire $J_c(\mathbf{B}, \theta) = \mathbf{J}_c$ et $n(\mathbf{B}, \theta) = \mathbf{n}$. Pour prendre en compte les dépendances en champ magnétique **B** et température $\theta$, il est nécessaire d'introduire des lois d'évolutions de ces paramètres, comme par exemple le modèle de Kim [24], et de combiner la résolution du problème magnétique à celle de l'équation de la chaleur. Dans ce cas, les para-
mètres $J_c$ et n sont calculés explicitement à chaque pas de temps en utilisant les solutions du pas de temps précédent [RI.6].

La linéarisation de la loi constitutive proposée est une approximation au premier ordre de la relation $\mathbf{E} = \rho(\mathbf{J})\mathbf{J}$. Nous la présentons sous la forme

$$ \mathbf{E}^k = \mathbf{E}^{k-1} + (\mathbf{J}^k - \mathbf{J}^{k-1}) \frac{\partial \mathbf{E}}{\partial \mathbf{J}} \quad (3.7) $$

avec k l'indice d'itération du processus de la linéarisation. Les champs $\mathbf{E}$ et $\mathbf{J}$ étant vectoriels, il est nécessaire de définir l'opérateur de dérivation vectoriel $\frac{\partial \mathbf{E}}{\partial \mathbf{J}}$. Nous avons proposé d'évaluer cet opérateur par un tenseur d'ordre 2 dont les termes sont les dérivées des composantes du champ électrique $\mathbf{E} = (E_1, E_2, E_3)$ par rapport aux composantes de la densité de courant $\mathbf{J} = (J_1, J_2, J_3)$ :

$$ \frac{\partial \mathbf{E}}{\partial \mathbf{J}} = \begin{pmatrix} \frac{\partial E_1}{\partial J_1} & \frac{\partial E_2}{\partial J_1} & \frac{\partial E_3}{\partial J_1} \\ \frac{\partial E_1}{\partial J_2} & \frac{\partial E_2}{\partial J_2} & \frac{\partial E_3}{\partial J_2} \\ \frac{\partial E_1}{\partial J_3} & \frac{\partial E_2}{\partial J_3} & \frac{\partial E_3}{\partial J_3} \end{pmatrix} \quad (3.8) $$

La résolution de la formulation en $\mathbf{H}$ en utilisant cette linéarisation a été implémentée dans une méthode des éléments finis. Elle permet une meilleure convergence notamment pour les cas où l'exposant de la loi puissance est grand. Ces résultats font l'objet d'un "benchmark" dans l'open source ONELAB développé par l'équipe du professeur Christophe Geuzaine de l'université de Liège en Belgique (https://gitlab.onelab.info/doc/models/wikis/Superconductors).

## 3.3 Méthodes avec flux pour le calcul de champs dans les supraconducteurs

L'approche des éléments finis est largement utilisée pour la modélisation des supraconducteurs. Elle présente néanmoins des points faibles tels que la consommation d'espace mémoire et le mauvais conditionnement des matrices à cause de la forte non-linéarité des équations. Il est alors important d'explorer d'autres méthodes notamment celles qui se sont avérées efficaces dans d'autres domaines, notamment en mécanique des fluides. Des travaux sur plusieurs méthodes ont été réalisés, comme par exemple, sur la méthode des différences finis et sur une méthode stochastique, qui ont été implémentées pour traiter des configurations 2D [11] [77]. Ma contribution au développement d'outils pour le calcul de champs dans ces matériaux a été de travailler sur l'utilisation des méthodes à "flux". J'ai notamment implémenté une méthode
des volumes finis et une méthode de Galerkin discontinue pour le calcul de la densité de courant induite dans les supraconducteurs.

### 3.3.1 Méthode des volumes finis

Mes travaux de thèse sur la modélisation des matériaux supraconducteurs se sont articulés autour du développement d'une méthode de volumes finis pour le calcul de la densité de courant induite dans le cas de configurations 2D [56]. Ces travaux ont montré que cette méthode est bien adaptée pour la résolution du problème de diffusion non-linéaire que l'on obtient lorsque le champ électrique n'a qu'une composante non nulle et satisfait à la condition $divE = 0$. Ce problème est donné par la formulation en E et s'écrit

$$ \left\{ \begin{array}{l} \mu_0 \frac{\partial J}{\partial t} - \Delta E = 0 \quad \text{in } \Omega \\ \nabla E \cdot \mathbf{v} = \nu_x \frac{\partial B_y}{\partial t} - \nu_y \frac{\partial B_x}{\partial t} \quad \text{on } \partial \Omega \end{array} \right. \qquad (3.9) $$

avec, E et J les inconnues scalaires du problème, $B = (B_x, B_y)$ le champ magnétique appliqué, $\mathbf{v} = (v_x, v_y)$ la normale extérieure au bord. Les configurations traitées sont celles de supraconducteurs massifs transportant du courant ou soumis à un champ magnétique transverse. Ces excitations sont prises en compte dans le modèle via une condition de Neumann non homogène qui est directement appliquée aux bords du matériau.

Soit $\Omega = \bigcup D$, un maillage du domaine de calcul en cellules élémentaires. Pour établir les équations algébriques à résoudre, on intègre l'équation différentielle sur chaque cellule $D_i$ en utilisant le théorème de la divergence, ce qui donne pour le système (3.9) :

$$ \mu_0 \int_{D_i} \frac{\partial J_i}{\partial t} dD_i - \int_{\Gamma} \nabla E_i \cdot \mathbf{v} d\Gamma_i \qquad (3.10) $$

avec $\Gamma_i = \partial D_i$, $J_i$ et $E_i$ respectivement la densité de courant et le champ électrique sur la cellule $D_i$ et $\mathbf{v}$ la normale extérieure à $\Gamma_i$. L'intégrale sur le bord $\Gamma_i$ est remplacée par une expression de flux numérique qui dépend de la solution sur les cellules voisines $D_j$ et prend en compte les propriétés de continuité :

$$ \int_{\Gamma_i} \nabla E_i \cdot \mathbf{v} d\Gamma_i \approx \sum_{D_j} |D_i \cap D_j| \Phi_{i,j} \qquad (3.11) $$

avec $|D_i \cap D_j|$ l'intersection entre les cellules $D_i$ et $D_j$ et $\Phi_{i,j}$ le flux numérique à construire. Dans la méthode des volumes finis construite, les expressions de flux numériques utilisées sont obtenues à partir des coefficients de la matrice de rigidité issue de la discrétisation de l'opérateur Laplacien par une méthode des éléments finis nodaux sur un maillage dual (Figure.3.1).
<figure><img src="image_58.png" /><figcaption>FIGURE 3.1 – (gauche) maillage éléments finis avec les noeuds au milieu des arêtes d’une triangulation de Ω, (droit) maillage volumes finis où les cellules sont construites autour des noeuds du maillage éléments finis Ω</figcaption></figure>

L'équation algébrique à résoudre pour une cellule Dᵢ à chaque pas de temps est la suivante :

$$ \mu_0 |D_i| \frac{J^p - J^{p-1}}{\delta t} - \sum_{D_j} S_{i,j} (E_i^p - E_j^p) = 0 \quad (3.12) $$

avec $J_i^p$ et $E_i^p$ respectivement la densité de courant et le champ électrique sur la cellule $D_i$ à l'instant $t_p$. La linéarisation proposée pour traiter la non-linéarité donnée par la relation (3.7) est introduite dans cette équation algébrique et une méthode de Newton Raphson est appliquée. On a ainsi pu traiter des cas où l'exposant de la loi puissance $n = 200$.

Ce travail a fait l'objet de publications dont celle-ci après **A 2-D Robust FE-FV Mixed Method to Handle Strong Nonlinearities in Superconductors** publiée dans le journal *IEEE Transactions on Magnetics*. On s'est appuyé sur ce modèle numérique pour mettre en évidence certains éléments sur la caractérisation de ces matériaux notamment celui de l'influence de la vitesse de variation du champ magnétique appliqué sur les pertes et l'aimantation des ces matériaux [RI.24].
# A 2-D Robust FE-FV Mixed Method to Handle Strong Nonlinearities in Superconductors

Abelin Kameni¹, Smail Mezani¹, Frédéric Sirois², Denis Netter¹, Jean Lévêque¹, and Bruno Douine¹

¹GREEN-UHP-INPL, Faculté des Sciences, 54506, BP. 239, Vandoeuvre-lès-Nancy, France

²Ecole Polytechnique de Montréal, Montréal, QC H3C 3A7, Canada

A robust numerical method based on 2-D mixed finite-elements-finite volumes (FE-FV) allows the solution of diffusion problems in superconducting (SC) materials. The proposed approach handles the strong nonlinearity of the E(J) constitutive power law of high-temperature superconductors (HTS). The method is tested for a SC cylinder submitted to a sinusoidal transport current or to a transverse sinusoidal external field. The current density distributions as well as the AC losses are computed. Comparisons to a FE analyses that use the magnetic field as state variable show the validity of the proposed approach. It can be seen that the proposed method is very stable even for large n-values for which the FE method does not converge.

**Index Terms**—Finite elements, finite volumes, superconducting materials.

## I. INTRODUCTION

THE CONSTITUTIVE power law is widely used to characterize high temperature superconductors. It is written as

$$J = J_c(E/E_c)^{1/n} \quad (1)$$

where *J* is the current density, *E* the electric field, *E*<sub>c</sub> the critical electric field, *J*<sub>c</sub> the critical current density, and *n* the power law exponent. *n* = 1 corresponds to a normal conductor and *n* = ∞ represents the critical state model suggested by Bean.

Finite-element methods are widely used to compute the magnetic field distribution in superconductors using different formulations [1], [2]. Unfortunately, the use of (1) leads to numerical oscillations and convergence problems for large n-values (typically, for *n* > 50).

In order to solve the nonlinear diffusion-convection-reaction equations in porous media, a coupled finite-elements-finite volumes (FE-FV) method has been successfully used [3]. In this paper, we propose to use a similar approach to solve the electric field diffusion equation in HTS materials. Applications to a superconducting cylinder under transport current or submitted to an external magnetic field are considered. It has been shown that the proposed approach is robust even for large n-values for which conventional approaches based on FE methods do not converge.

## II. PRINCIPLE OF THE METHOD

We set $u = E/E_c$, $\beta(u) = J/J_c$, and $c = \mu_0 J_c / E_c$.

The superconductor has a vacuum magnetic permeability $\mu_0$.

In 2-D problems, *J* and *E* are scalar quantities having only one component in the *z* direction. The SC domain is noted $\Omega$ and its external boundary $\Gamma$. Using Maxwell equations and the constitutive law (1), the problem to solve is

$$\left\{ \begin{array}{ll} c \frac{\partial^3 u}{\partial t^3} - \Delta u = 0 & \text{in } \Omega \\ \nabla \cdot \vec{u} \cdot \vec{\nu} = C_b(\nu) & \text{on } \Gamma \end{array} \right. \quad (2)$$

Manuscript received December 21, 2009; revised February 11, 2010; accepted February 16, 2010. Current version published July 21, 2010. Corresponding author: S. Mezani (e-mail: smail.mezani@green.uhp-nancy.fr).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TMAG.2010.2044025

The boundary condition on $\Gamma$ results from the variation of the magnetic field according to Faraday's law

$$C_b(\nu) = E_c^{-1} (\partial B_y / \partial t_r - \partial B_r / \partial t_s) \cdot \vec{\nu} \quad (3)$$

$B_r, B_y$ are, respectively, the *x*, *y* components of the flux density and $\vec{\nu}$ is the outward normal vector on $\Gamma$.

The discretization methods (FEM or FVM) are applied to the weak formulation of (2) in the time interval $[t_p, t_{p+1}]$

$$\int_{\Omega} c \frac{\beta(u^{p+1}) - \beta(u^p)}{\Delta t_p} \varphi d\Omega - \int_{\Omega} \Delta u^{p+1} \varphi d\Omega = 0 \quad (4)$$

where $\varphi$ is a test function, $u^p$ is the solution at instant $t_p$, and $\Delta t_p$ the time step.

The method exposed here is based on the FE-FV coupling principle. It exploits the duality of the FE and FV meshes as well as the footbridge operators between FE and FV approximation spaces, and combines FE and FV discretizations terms of (4).

### A. FE-FV Coupling Principle

The FE-FV coupling is widely used in fluid mechanics to solve convection-diffusion type equations. The computation procedure consists in successively solving the diffusion equation on the FE mesh and the convection equation on the FV cell. In so doing, it is necessary to calculate the projection of the FE solution on the FV mesh and vice versa.

Starting from a solution $u^p$ at time $t_p$, the solution $u^{p+1}$ at time $t_{p+1}$ using the FF-FV coupling scheme is determined as follows:

i) calculate the solution $u_{p,1}$ at $t_p$ of the FV approach;

ii) determine its projection $u_{p,2}$ on the FE mesh;

iii) use $u_{p,2}$ to compute the solution $u_{p,3}$ at $t_{p+1}$ of the FE formulation;

iv) the needed solution $u_{p+1}$ is the projection of $u_{p,3}$ on the FV mesh.

The problem treated here is a convection-diffusion equation. The objective is to combine the advantages of each method (FE and FV) in a single numerical scheme to solve the problem given in (2). In so doing, we have to ensure that:

- the two approaches give the same solution;

- the diffusion on a FE node is equivalent to the convection on the corresponding FV cell.
<figure><img src="image_59.png" /><figcaption>Fig. 1. Association of each node of the FE mesh to the cells of the FV mesh: we link the centers of gravity of two contiguous triangles through the vertex of their common edge. (a) The nodes Ni correspond to the middle of the edges. (b) The nodes Ni correspond to the center of the cells $D_i$.</figcaption></figure>

## B. Footbridge Operators

Let $u^{\text{FE}}$ and $u^{\text{FV}}$ be the respective solutions of the FE and FV approaches. The footbridge operators $\pi_{\text{E}}(\pi_{\text{E}}(u^{\text{FV}})) = u^{\text{FE}}$ and $\pi_{\text{V}}(\pi_{\text{V}}(u^{\text{FE}})) = u^{\text{FV}}$ between the FE and the FV approximation spaces are well described in [4]. They are defined as

$$ \int_{\Omega} |\pi_{\text{E}}(u) - u|^2 d\Omega = \min_{u \in \Gamma_{\text{FE}}} \int_{\Omega} |w - u|^2 d\Omega \quad (5) $$

$$ \int_{\Omega} \pi_{\text{V}}(u) \varphi d\Omega = \int_{\Omega} u \varphi d\Omega. \quad (6) $$

Let $\mathfrak{S}$ a regular triangulation of $\Omega$ and $\vartheta$ a FV partition of $\Omega$. When $\vartheta$ is obtained from $\mathfrak{S}$, the duality between the FE and FV meshes is described by the identification of each FE node to each FV cell, Fig. 1. The basis functions of the FV approximation space are obtained from the indicatrix functions of each cell of the FV mesh. We remark that each first order FE base functions $\varphi_i$ constructed on the nodes $N_i$ of the FE triangulation can be prolonged on the corresponding FV cell.

In the case of the Fig. 1, it has been proved in a previous paper [4] that the footbridge operators ensure the equality $u^{\text{FE}} = u^{\text{FV}}$.

## C. Numerical Scheme

The numerical discretization of the Laplacian operator is well established using the FE method whereas it is not adapted to the strong nonlinearity of the $\beta$ term. Since $\beta$ is continuous, its mean value is well defined so we use the FV method to deal with the nonlinearity

$$ \int_D \frac{\beta(u^{p+1}) - \beta(u^p)}{\Delta t_p} dD = \frac{\beta^{p+1} - \beta^p}{\Delta t_p} D \quad (7) $$

where $D$ is the surface of a FV cell and $\beta^p$ is the mean value of $\beta(u^p)$.

Hence, using the footbridge operators, we combine FE-FV methods to take advantage of both of them in the same numerical scheme.

In the FV method, the flux $\Phi_{i,j}$ between two adjacent cells $i$ and $j$ must be continuous. The mixed numerical scheme is justified if the diffusion on a FE node is equal to the convection in the FV cell. This leads to

$$ \underbrace{\int_{D_i} \operatorname{div}(\vec{\nabla} u_i) dD_i}_{\Gamma_{\text{FV}}} = \underbrace{\int_{\Omega} \vec{\nabla} u_i \cdot \vec{\nabla} \varphi_i d\Omega}_{\Gamma_{\text{FE}}}. \quad (8) $$

To ensure the equality in (8), the $\Phi_{i,j}$ are constructed from the stiffness matrix coefficient $S_{i,j}$ of the FE discretization as follows:

$$ \sum_{j \neq i} |D_i \cap D_j| \Phi_{i,j} = - \sum_{j} S_{i,j} u_j \quad (9) $$

where $s_{i,j} = \sum_{K \in \mathfrak{S}} \int_K \nabla \varphi_i \cdot \nabla \varphi_j dK$.

If we define $S_{i,j} = - \sum_{j \neq i} S_{i,j}$, the local conservation property of the FV formulation is ensured.

The discrete problem is formed using (7) and (9). For the internal FV cells, the discrete problem is written as

$$ c_i \frac{\beta_i^{p+1} - \beta_i^p}{\Delta t_p} D_i = \sum_j S_{i,j} u_j^{p+1}. \quad (10) $$

and for the cells on the boundary $\Gamma$, we write

$$ c_i \frac{\beta_i^{p+1} - \beta_i^p}{\Delta t_p} D_i = \sum_j S_{i,j} u_j^{p+1} + |\Gamma_i| C_b^{p+1}. \quad (11) $$

The existence of the discrete solution is proved by the saddle-point theorem. The numerical computation uses the unknowns $u_i^{p+1} = \beta(u_i^{p+1})$. Since $\beta^{-1}$ is continuous and differentiable, we can use a Newton-Raphson algorithm

$$ u_i^{p,k} = \beta^{-1} (u_i^{p,k-1}) + (\beta^{-1})' (u_i^{p,k-1}) \times (u_i^{p,k} - u_i^{p,k-1}). \quad (12) $$

## III. APPLICATION EXAMPLES

We consider a superconducting cylinder having a radius $R = 1.5$ mm and characterized by $J_c = 14.15$ A/mm$^2$, $E_c = 10^{-4}$ V/m. Two values of $n$ are considered; $n = 25$ and $n = 200$. The first value is typical for high temperature superconductors (BiSCCO or YBCO materials) and the second corresponds to low temperature superconductors (NbTi for example). For this last case, the Bean model can be used to validate the proposed method.

### A. Cylinder Under Self Field (Transport Current)

The cylinder carries a sinusoidal transport current $i(t) = I_{\text{max}} \sin(2\pi f t)$, $I_{\text{max}} = 40.8$ A. Two frequency values are considered, $f = 0.5$ Hz and $f = 1$ kHz.

$B_r, B_y$ used in the boundary condition (3) are obtained by Ampere's theorem.

Fig. 2 presents the current density distribution on the $Ox$ axis at $t = T/4$ (T is the period) for the two considered frequencies. A comparison to the results issued from a purely FE formulation [2] where the magnetic field is used as state variable is also given. A good agreement is observed.

We have also compared in Fig. 3, the instantaneous self field losses obtained by the FE-FV method and those issued from the FE method. Once again, satisfactory results are obtained. These results validate the proposed FE-FV approach for weak $n$ values.
<figure><img src="image_60.png" /><figcaption>Fig. 2. Comparison of current density profiles for $n = 25$. (a) $f = 0.5$ Hz at $t = T/4 = 0.5$ s. (b) $f = 1$ kHz at $t = T/4 = 0.25$ ms.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>Fig. 3. Comparison of instantaneous self field losses for $n = 25$, $f = 0.5$ Hz.</figcaption></figure>

In the case where $n = 200$, the FE method does not converge. The average AC losses obtained by our method are 0.024 W/mm. The Bean formula [6] gives 0.028 W/mm which also validate the FE-FV method in the case of large $n$ values.

<figure><img src="image_62.png" /><figcaption>Fig. 4. Current density distribution $J/J_c$ at $t = T/4$ for $n = 25$ ($B_{max} = 15$ mT). (a) $f = 1$ kHz. (b) $f = 0.5$ Hz.</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>Fig. 5. Current density distribution $J/J_c$ at $t = T/4$ for $n = 200$ ($B_{max} = 15$ mT). (a) $f = 1$ kHz. (b) $f = 0.5$ Hz.</figcaption></figure>

## B. Cylinder Submitted to an External Magnetic Field

The cylinder is submitted to a sinusoidal transverse magnetic field in the $x$ direction, $b_r(t) = B_{max} \sin m(2\pi f t)$. Two frequency values are considered, $f = 0.5$ Hz and $f = 1$ kHz.

The current density distribution in the cylinder is presented in Fig. 4 for $n = 25$ and Fig. 5 for $n = 200$. The results at $t = T/4$ (when the applied field reaches its maximum value) and $n = 25$ show that the penetration of the current density at $f = 0.5$ Hz is more important than the penetration at $f = 1$ kHz. However, the current density reaches a value $J/J_c = 0.83$ for $f = 0.5$ Hz and $J/J_c = 1.3$ for $f = 1$ kHz. This shows the influence of the speed variation (i.e., the frequency) of the applied field.

For $n = 200$, it can be seen that the penetration of the current density is almost the same for both frequencies. This is in accordance with the Bean's model for which the speed variation of the applied field does not play any role when $n \to \infty$.

Fig. 6 presents the magnetization versus applied field. In the simulation, $B_{max}$ is equal to 50 mT, a value for which the cylinder is in a complete penetration state.

For $n = 25$, Fig. 6(a), the maximum value of the magnetization is higher for $f = 1$ kHz than the value obtained at $f = 0.5$ Hz. Again, this shows the influence of the speed variation (i.e., the frequency) of the applied field.

For $n = 200$, the influence of the frequency is practically inexistent. The maximum value of the magnetization given by the Bean model can be computed as $M_{max} = J_c \times R/3$, which leads to a value $M_{max} = 7.07$ A/mm for the considered example. Using the proposed method, Fig. 6(b), we obtain $M_{max} = 6.95$ A/mm. Again, this validates our method in the case of large $n$ values.
<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /><figcaption>Fig. 6. Magnetization against applied field ($B_{max}$ = 50 mT). (a) $n$ = 25. (b) $n$ = 200.</figcaption></figure>

## IV. CONCLUSION

In this paper, we proposed a mixed finite elements-finite volumes numerical scheme to solve the 2-D scalar diffusion of the electric field in superconducting materials. It has been shown that the proposed approach can successfully handle the strong nonlinearity of the E-J power law. Comparisons to purely FE analyses show the validity of our method in the case of weak n-values. In the case of large n-values for which the FE method does not converge, it has been shown via the comparisons to Bean's model that the proposed approach is very robust to treat the strong nonlinearities.

## REFERENCES

[1] E. Vinot, G. Meunier, and P. Tixador, "Different formulations to model superconductors," *IEEE Trans. Magn.*, vol. 36, no. 4, pp. 1226-1229, Jul. 2000.

[2] R. Brambilla, F. Grilli, and L. Martini, "Development of an edge-element model for AC losses computation of high-temperature superconductors," *Supercond. Sci. Technol.*, vol. 19, pp. 1-9, 2006.

[3] R. Eymard, D. Hilhorst, and M. Vohralík, "A combine finite volume-nonconforming/mixed-hybrid finite element scheme for degenerate parabolic problems," *Numer. Math.*, vol. 105, pp. 73-131, 2006.

[4] F. Pascal and J. M. Ghidaglia, "Footbridge between finite volumes and finite elements with applications to CFD," *Int. J. Numer. Methods Fluids*, vol. 37, no. 8, pp. 951-986, 2001.

[5] K. Kajikawwa, Y. Fukuda, H. Moriyama, T. Hayashi, M. Iwakuma, and K. Funaki, "Frequency dependence of AC magnetization in stacked Bi-2223 Ag-sheathed tapes," *Physica C*, vol. 382, pp. 122-126, 2002.

[6] W. J. Carr, Jr., *AC Loss and Macroscopic Theory of Superconductors*. London, U.K.: Taylor & Francis, 2001.
3.3.2 Méthode de Galerkin discontinue

La méthode de Galerkin discontinue est de plus en plus employée dans des travaux de mo-
délisation en électromagnétisme. Son utilisation pour la modélisation des matériaux supracon-
ducteurs a été introduite grâce à nos travaux [RI.6, RI.9, RI.17, RI.23] et de nouveaux déve-
loppements ont récemment été proposés [26]. Cette méthode étant bien adaptée à la mise en
place du calcul parallèle et à l'utilisation des éléments d'ordre élevé, offre plusieurs possibilités
d'amélioration des modèles numériques autour de ces matériaux. Ainsi, les coûts de simulation
peuvent être réduits grâce au calcul parallèle et la précision des résultats peut être améliorée
grâce à l'utilisation des éléments d'orde élevé. Nous avons mis en place cette méthode pour la
résolution des formulations en E et H données respectivement par les relations (3.3) et (3.2).

Résolution de la formulation électrique

Dans mes travaux de thèse j'ai implémenté une méthode des volumes finis pour la résolution
de la formulation en **E** pour des configurations 2D où il n'y a pas de champs démagnétisants
et où l'inconnue du problème n'a qu'une composante non-nulle. Dans le cas plus général, l'in-
connue est vectorielle et les effets démagnétisants doivent être pris en compte.

Nous avons proposé d'étendre les travaux de thèse, à des cas où le champ électrique **E** a ses trois composantes non-nulles. Pour cela, nous avons supposé que le problème vectoriel satisfait par **E** (3.3) est la superposition de trois équations scalaires satisfaites par les composantes (E₁, E₂, E₃). Pour chacune des inconnues Eᵢ, nous avons construit une condition au bord de type Neumann non homogène en utilisant un développement explicite de la relation de Maxwell-Faraday [RI.23]. Le système pour chacune de composante est de la forme :

$$
\left\{
\begin{array}{l}
\mu_0 \frac{\partial J_i}{\partial t} - \Delta E_i = S_i \quad \text{in } \Omega \\
\nabla \cdot E_i \cdot \mathbf{v} = f \left( \frac{\partial \mathbf{B}}{\partial t}, \mathbf{E} \right) \quad \text{on } \partial \Omega
\end{array}
\right.
\qquad (3.13)
$$

avec i ∈ {1, 2, 3}, Sᵢ est la ième composante de grad(divE), et f une fonction qui dépend de
la variation du champ magnétique total et des solutions Eᵢ calculées au pas de temps précé-
dent. Cette fonction représente le flux au bord qui est imposé directement sur le matériau. Pour
éviter d'avoir à ajouter un volume d'air autour du matériau pour prendre en compte les effets
démagnétisants, nous avons proposé d'utiliser la relation qui donne le champ démagnétisant
en fonction de l'aimantation et les coefficients démagnétisants. Ce champ démagnétisant est
calculé à chaque itération et est introduit dans les conditions aux bords. Cette approche pour
prendre en compte les effets démagnétisants est uniquement utilisable pour des géométries dont
ont connait les coefficients démagnétisants [RI.17].
La discrétisation par la méthode Galerkin discontinue du système (3.13) s'appuie sur la formu-
lation faible établie sur chaque élément K du maillage $\Omega = \bigcup K$ :

$$
\int_K \mu_0 \frac{\partial J_i}{\partial t} \varphi dK - \int_K \nabla E_i \cdot \nabla \varphi dK - \int_{\Gamma_K} \nabla E_i \cdot \nu \varphi d\Gamma_K = \int_K S_i \varphi dK \quad (3.14)
$$

Elle consiste à appliquer une approche éléments finis à valeurs nodales pour discrétiser la for-
mulation faible sur chaque élément K (3.14) en utilisant des expressions de flux numériques
pour remplacer les intégrales de bords comme dans une méthode des volumes finis. Les expres-
sions de flux numériques adoptées pour ce système, contiennent des termes de pénalités per-
mettant de satisfaire aux propriétés de continuité. Dans notre travail, nous avons implémenté les
expressions de flux numériques issues, de la méthode de pénalité intérieure symétrique (Sym-
metric Interior Penalty method - SIP) et de la méthode de pénalité intérieure non-symétrique
(Non-symmetric Interior Penalty method - NIP), dont les expressions sont :

$$
\nabla E_i \cdot \nu \rightarrow \begin{cases} \{\nabla E_i\} \cdot [[\varphi]]\nu + \alpha [[E_i]][[\varphi]] & \text{(NIP)} \\ [[\nabla E_i]] \cdot \{\varphi\}\nu + \{\nabla E_i\} \cdot [[\varphi]]\nu + \alpha [[E_i]][[\varphi]] & \text{(SIP)} \end{cases} \quad (3.15)
$$

où $\alpha$ est le coefficient de pénalité, $\{\bullet\}$ et $[[\bullet]]$ désignent respectivement la fonction moyenne et la fonction saut à l'interface.

Ce travail à fait l'objet de deux publications **Discontinuous Galerkin Method for Computing Induced Fields in Superconducting Materials** et **A 3-D Semi-Implicit Method for Computing the Current Density in Bulk Superconductors** publiées dans *IEEE Transactions On Magnetics*.
# Discontinuous Galerkin Method for Computing Induced Fields in Superconducting Materials

A. Kameni<fnref n="1" />,<fnref n="2" />, J. Lambrechts<fnref n="3" />, J.-F. Remacle<fnref n="3" />, S. Mezani⁴, F. Bouillault<fnref n="2" />, and C. Geuzaine<fnref n="1" />

A discontinuous Galerkin method is proposed for computing the current density in superconductors characterized by a constitutive power law between the current density and the electric field. The method is formulated to solve the nonlinear diffusion problem satisfied by the electric field, both in the scalar and 2-D vectorial case. Application examples are given for a superconducting cylinder subjected to an external magnetic field. Results are compared to those given by the mixed finite-element/finite-volume method and those obtained using a standard finite-element software. Efficiency and robustness of the approach are illustrated on an example where the exponent in the power law is spatially dependent.

**Index Terms**—Discontinuous Galerkin method, finite-element method, finite-volume method, nonlinear diffusion, superconductors.

## I. INTRODUCTION

THE constitutive power laws $E - J$ and $J - E$ are widely used to characterize high temperature superconductors. They are written as

$$ \frac{\vec{J}}{J_c} = \left\| \frac{\vec{E}}{E_c} \right\|^{\frac{1}{n}-1} \frac{\vec{E}}{E_c} \quad (1) $$

$$ \frac{\vec{E}}{E_c} = \left\| \frac{\vec{J}}{J_c} \right\|^{n-1} \frac{\vec{J}}{J_c} \quad (2) $$

where $\vec{J}$ is the current density, $\vec{E}$ the electric field, $E_c$ the critical electric field, $J_c$ the critical current density, and $n$ the power law exponent. The case $n = 1$ corresponds to a normal conductor, while $n = +\infty$ represents the critical state model suggested by Bean [1]. Several numerical methods have been proposed to solve nonlinear diffusion problems resulting from Maxwell's equations [2]–[5]. Their results are satisfying when $n$ is uniform and sufficiently small. Few of them are suited when $n$ is large, and models where $n$ locally varies are uncommon.

In this paper, we present a Discontinuous Galerkin (DG) method for computing induced fields in superconductors. We work on solving the nonlinear diffusion problem in terms of the electric field in order to determine the current density when $n$ is large or locally varies. DG methods are well suited to treat discontinuous forms and use a piecewise high-order polynomials basis for reducing spurious oscillations. In addition, they are naturally well adapted for parallel computing.

Manuscript received July 07, 2011; revised October 05, 2011; accepted October 21, 2011. Date of current version January 25, 2012. Corresponding author: C. Geuzaine (e-mail: cgeuzaine@ulg.ac.be).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TMAG.2011.2173664

## II. THE DIFFERENTIAL SYSTEM

In a two-dimensional setting where the magnetic induction depends only on two space variables ($\vec{B} = (B_r, B_y)$), the electric field $\vec{E}$ and current density $\vec{J}$ have a single nonzero component and can thus be treated as scalar fields $E$ and $J$. We set $u = E/E_c$, $\beta(u) = u^{1/n} = J/J_c$, and $c = \mu_0 J_c/E_c$. The superconductor sample has a vacuum magnetic permeability $\mu_0$. Denoting the superconductor domain by $\Omega$ and its border by $\partial\Omega$, Maxwell's equations and the constitutive law (1) lead to the following nonlinear diffusion problem:

$$ (S) \left\{ \begin{array}{l} \frac{\partial \beta(u)}{\partial t} - \frac{1}{c} \Delta u = 0 \quad \text{in } \Omega \\ \vec{\nabla} \cdot \vec{J} \cdot \vec{u} = C_b(t) \quad \text{on } \partial \Omega. \end{array} \right. \quad (3) $$

The system is established with a zero initial condition and the boundary condition on $\partial\Omega$ results from Faraday's law:

$$ C_b(t) = E_c^{-1} \left( \frac{\partial B_y}{\partial t} - \frac{\partial B_r}{\partial t} \right) \cdot \vec{v} \quad (4) $$

where $\vec{v}$ is the outward normal vector.

## III. DISCONTINUOUS GALERKIN METHOD

Let us consider a triangulation $\mathcal{T}_h = \bigcup K$ of the domain $\Omega$. The DG approach combines discretization tools of finite-element (FE) and finite-volume (FV) methods. It consists in solving on each $K$ the weak formulation of the system (3):

$$ c \int_K \frac{\partial \beta(u)}{\partial t} \varphi dK - \int_K \vec{\nabla}_u \cdot \vec{\nabla} \varphi dK - \int_{\Gamma} \vec{\nabla}_u \cdot \vec{v} \varphi d\Gamma = 0 \quad (5) $$

where $\varphi$ is a test function and $\Gamma = \partial K$ is an interface between two elements of $\mathcal{T}_h$ or a part of $\partial\Omega$.

### A. FE Discretization on Each Element

On each triangle $K \in \mathcal{T}_h$ a finite-element approximation space is defined. Its basis functions are Lagrange polynomials

<footnote n="1">¹Université de Liège, Department of Electrical Engineering and Computer Science, Montefiore Institute, Liège 91192, Belgium</footnote>

<footnote n="2">²Laboratoire de Génie Electrique de Paris, Gif-sur-Yvette 91192, France</footnote>

<footnote n="3">³Université Catholique de Louvain, Center for Systems Engineering and Applied Mechanics, Louvain-la-Neuve B-1348, Belgium</footnote>

<footnote n="1">¹Groupe de Recherche en Electronique et Electrotechnique de Nancy, Vandoeuvre-la-Nancy 54506, France</footnote>
of degree $p$. The number of nodes on $K$ is given by $(1/2)(p+1)(p+2)$. The discrete solution is written as

$$u^K = \sum_j u_j^K \varphi_j^K. \qquad (6)$$

with $u_j^K$ its value at node $j$. Since $\beta$ is a Liotschitz function, we assume that

$$\beta(u^K) = \sum_j \beta(u_j^K) \varphi_j^K. \qquad (7)$$

The elementary mass matrix $M^K$ is determined from the $L^2$ scalar product:

$$M_{ij}^K = \int_K \varphi_i^K \varphi_j^K dK \qquad (8)$$

and the global mass matrix $M$ is block diagonal.

## B. Flux Term on the Interface Between Two Elements

As in an FV method, the interface term $\int_{\Gamma} \vec{\nabla} u \cdot \vec{\nu} \varphi d\Gamma$ is treated by a numerical flux $F$, which verifies $F^{K,L} = -F^{L,K}$, where $K$ and $L$ are neighboring elements. Its construction needs the following functions at the interface $K \cap L$: the mean value

$$\{u\} = \frac{u^L + u^K}{2} \qquad (9)$$

and the jump

$$[|u|] = u^L - u^K. \qquad (10)$$

Many expressions of $F$ have been proposed in the case of the Laplace operator. We choose the expression of $F$ based on the Non-symmetric Interior Penalty method (NIP). The NIP method consists of introducing a penalty term $\int_{\Gamma} \theta[|u|][|\varphi|]d\Gamma$, in order to guarantee continuity of $u$ and $\vec{\nabla} u$ at the interface [6]. The numerical flux is given by

$$F \rightarrow \int_{\Gamma} \{\vec{\nabla} u\} \cdot [|\varphi|] \vec{\nu} d\Gamma + \int_{\Gamma} \theta [|u|] [|\varphi|] d\Gamma \quad (11)$$

where $\theta$ is a positive parameter.

## C. The Discrete Problem

Rules for evaluating the different terms of the weak formulation (5) exploit properties of mesh parametrization also called “mapping”. The mapping is based on a bijective function $\Psi$ such as $\Psi(x, y) = (\xi, \eta)$, which allows to transform the physical space $(x, y)$ into a (reference) parametric space $(\xi, \eta)$ [7] (see Fig. 1). The basis functions in the parametric space are the linear combinations of $\xi^i \eta^j$, where $\alpha + \beta \le p$. The terms of the weak formulation (5) are evaluated in the parametric space, where derivation and integration operations are more convenient, and then mapped in the physical space.

After time discretization a discrete problem is obtained on each triangle $K$:

$$M^K \frac{\beta(u_{k+1}^K) - \beta(u_k^K)}{\delta t} = f^K(u_{k+1}^K) \qquad (12)$$

<figure><img src="image_66.png" /><figcaption>Fig. 1. Mapping of a triangle, from the (reference) parametric space to the physical space.</figcaption></figure>

where $\delta t$ is the time step, $u_k^K$ is the solution at instant $t_k$, and $f^K(u_{k+1}^K)$ represents the discretization of the Laplace operator with numerical flux given in (11).

Unfortunately $\beta'(0) \to +\infty$ and we cannot directly use a Newton iterative method to solve (12). We thus use the change of variables: $v = \beta(u)$ and $g^K(v) = f^K(\beta^{-1}(v))$, so that the discrete problem (12) becomes

$$M^K \frac{v_{k+1}^K - v_k^K}{\delta t} = g^K(v_{k+1}^K) \qquad (13)$$

which can be solved using a Newton iterative method since $g^K$ is continuous and differentiable.

Considering the contributions of all the elements in the mesh (dropping the $K$ superscript), we introduce the global unknown vector $d\nu = v_{k+1} - v_k$, which will be evaluated at each iteration ($l$) of the Newton method, and we define $h$ such as

$$h(d\nu) = \delta t g(v_k + d\nu) - M d\nu = 0. \qquad (14)$$

The iterative solution $d\nu^{(l+1)}$ is obtained using the following relation:

$$\frac{\partial h(d\nu^{(l)})}{\partial d\nu} (d\nu^{(l+1)} - d\nu^{(l)}) = -h(d\nu^{(l)}). \qquad (15)$$

Thanks to relation (14), we have

$$\frac{\partial h(d\nu^{(l)})}{\partial d\nu} = \delta t g'(v_k + d\nu^{(l)}) - M \qquad (16)$$

so that equation (15) can be written as the following linear system:

$$\begin{aligned} & [\delta t g'(v_k + d\nu^{(l)}) - M] d\nu^{(l+1)} \\ & = [\delta t g'(v_k + d\nu^{(l)})] d\nu^{(l)} - \delta t g(v_k + d\nu^{(l)}) \end{aligned} \qquad (17)$$

where $g'(v_k + d\nu^{(l)}) = (\beta^{-1})'(v_k + d\nu^{(l)}) f'(\beta^{-1}(v_k + d\nu^{(l)}))$.

The Newton iterations are started with $d\nu^{(0)} = 0$ and stopped when $||h(d\nu^{(l)})|| < 1e-10$.

## IV. EXTENSION TO A SEMI-IMPLICIT SCHEME FOR THE 2-D VECTORIAL CASE

Let us consider an infinitely long superconductor along the $z$ axis. In a two-dimensional setting where the magnetic induction is axial ($\vec{B} = (0, 0, B_z)$), the electric field and current density have two nonzero components. The diffusion problem satisfied by the electric field becomes vectorial:

$$\frac{\partial \vec{J}}{\partial t} - \frac{1}{c} \vec{\nabla} \vec{E} = 0. \qquad (18)$$
Our idea for getting the solution of this vectorial problem is to solve independently the diffusion equations of each components $E_r, E_y$ of the electric field.

Setting $u_1 = E_r/E_c, u_2 = E_y/E_c$ and

$$ \beta_1(u_1, u_2) = (u_1^2 + u_2^2)^{\frac{1-\nu}{2\nu}} \quad u_1 = J_r / J_c = v_1 \quad (19) $$

$$ \beta_2(u_1, u_2) = (u_1^2 + u_2^2)^{\frac{1-\nu}{2\nu}} \quad u_2 = J_y / J_c = v_2 \quad (20) $$

the diffusion equations for $E_r$ and $E_y$ are written

$$ (S) \begin{cases} \frac{\partial \beta_1(u_1, u_2)}{\partial t} - \epsilon^{-1} \Delta u_1 = 0 & \text{in } \Omega \\ \frac{\partial \beta_2(u_1, u_2)}{\partial t} - \epsilon^{-1} \Delta u_2 = 0 & \text{in } \Omega \\ \vec{\nabla} u_1 \cdot \vec{\nu} = C_{b1}(t) & \text{on } \partial \Omega \\ \vec{\nabla} u_2 \cdot \vec{\nu} = C_{b2}(t) & \text{on } \partial \Omega \end{cases} \quad (21) $$

where $C_{b1,2}$ are built from Faraday's law and the hypothesis that $\vec{\nabla} \cdot \vec{E} = 0$ on $\partial \Omega$.

We propose a semi-implicit approach in order to handle the nonlinearity of the problem. At each time step, the explicit form $u_2^k$ is used to determine the solution $u_1^{k+1}$ (and vice versa). This assumption leads to the following expression of the nonlinearities:

$$ \beta_1(u_1^{k+1}, u_2^k) = \left[ (u_1^{k+1})^2 + (u_2^k)^2 \right]^{\frac{1-\nu}{2\nu}} u_1^{k+1} \quad (22) $$

$$ \beta_2(u_1^k, u_2^{k+1}) = \left[ (u_1^k)^2 + (u_2^{k+1})^2 \right]^{\frac{1-\nu}{2\nu}} u_2^{k+1} \quad (23) $$

The Neumann boundary conditions are constructed thanks to this semi-implicit formulation. In order to obtain the fluxes on the boundary, we consider that

$$ \vec{\nabla} \cdot \vec{E} = 0 \Rightarrow \begin{cases} \frac{\partial u_1^{k+1}}{\partial r} = - \frac{\partial u_2^k}{\partial y}, & \text{unknown is } u_1 \\ \frac{\partial u_1^{k+1}}{\partial y} = - \frac{\partial u_1^k}{\partial r}, & \text{unknown is } u_2 \end{cases} \quad (24) $$

$$ \vec{\nabla} \times \vec{E} = - \partial_t \vec{B} \Rightarrow \begin{cases} \frac{\partial u_1^{k+1}}{\partial y} = - \frac{\partial B}{\partial t} - \frac{\partial u_2^k}{\partial r}, & \text{unknown is } u_1 \\ \frac{\partial u_2^{k+1}}{\partial r} = \frac{\partial B}{\partial t} - \frac{\partial u_1^k}{\partial y}, & \text{unknown is } u_2 \end{cases} \quad (25) $$

where $C_{b1,2}$, given in (26) and (27), are deduced from the scalar product with the outward normal vector $\vec{\nu}$

$$ C_{b1}^{k+1} = - \frac{\partial u_2^k}{\partial y} \nu_r - \left( \frac{\partial B}{\partial t} + \frac{\partial u_2^k}{\partial x} \right) \nu_y \quad (26) $$

$$ C_{b2}^{k+1} = \left( \frac{\partial B}{\partial t} - \frac{\partial u_1^k}{\partial y} \right) \nu_r - \frac{\partial u_1^k}{\partial x} \nu_y \quad (27) $$

The discrete system to solve becomes

$$ \begin{cases} M K \cdot \beta_1(u_1^{k+1}, u_2^k) - \beta_1(u_1^k, u_2^k) = f_1(u_1^{k+1}, u_2^k) \\ M K \cdot \beta_2(u_1^k, u_2^k)^{\frac{\epsilon t}{\epsilon}} - \beta_2(u_1^k, u_2^k) = f_2(u_1^k, u_2^{k+1}) \end{cases} \quad (28) $$

where $f_1$ and $f_2$ are issued from the spatial discretization of the Laplace operator. Unfortunately the inverses of these semi-implicit functions $\beta_1$ and $\beta_2$ are nontrivial.

In accordance with the $E-J$ power law (2) we suppose that

$$ u_1^{k+1} = \beta_1^{-1}(v_1^{k+1}, v_2^k) = \left[ (v_1^{k+1})^2 + (v_2^k)^2 \right]^{\frac{1-\nu}{2}} v_1^{k+1} \quad (29) $$

<figure><img src="image_67.png" /><figcaption>Fig. 2. Current density distribution at $t = 0.5$ s with $n = 200$: (left) DG method with $p = 1$; (right) mixed FE-FV method.</figcaption></figure>

$$ u_2^{k+1} = \beta_2^{-1}(v_1^k, v_2^{k+1}) = \left[ (v_1^k)^2 + (v_2^{k+1})^2 \right]^{\frac{\nu-1}{2}} v_2^{k+1} \quad (30) $$

These inverse functions are continuous and differentiable. The previously described Newton iterative method is applied for computing $v_1^{k+1}$ and $v_2^{k+1}$.

## V. NUMERICAL RESULTS

We consider a superconducting cylinder of radius $R = 1.5$ mm, characterized by $J_c = 14.15$ A/mm$^2$ and $E_c = 10^{-1}$ V/m. The external magnetic flux density field is sinusoidal of amplitude $B_0$, frequency $f$ and period $T = 1/f$.

### A. Scalar Case

When the cylinder is subjected to an external transverse magnetic field in the $x$ direction the diffusion problem of the electric field is scalar and given by (3). We report numerical results obtained with $f = 0.5$ Hz, and $B_0 = 15$ mT.

For $n = 200$, the current density distribution is plotted at $t = T/4$. Notice that we did not get convergence using the $H$-formulation [2] implemented in the Comsol software. Fig. 2 presents a comparison of our result to that given by the mixed FE-FV method [3]. A good agreement is observed.

### B. The 2-D Vectorial Case

Let us now assume that the external magnetic field is applied in the $z$ direction. The diffusion problem vectorial and given by (21). In this example we consider $f = 50$ Hz and $T = 20$ ms. Based on the Bean model, the field required for full penetration will occur at $t = T/4$, where the applied field reaches its maximum value [8].

For $n = 200$, we present the distributions of the components $J_x/J_c, J_y/J_c$ of the norm $||\vec{J}||/J_c$ of the current density at different times. Fig. 3 shows the components in partial penetration at $t = T/8$. They have the expected symmetries because of similarities of their diffusion equations. The penetration form is different compared to the one obtained with a transverse magnetic field. Fig. 4 shows the components in full penetration at $t = T/4$. Symmetries of their distributions are unchanged. The regions where we get iso-values are in a pie chart form.

The modulus $||\vec{J}||/J_c = \sqrt{J_x^2/J_c^2 + J_y^2/J_c^2}$ is plotted in Fig. 5. In partial penetration at $t = T/8$, the repartition of the current density is orthoradial. The full penetration of the superconducting cylinder is reached at $t = T/4$ as in the Bean model.
<figure><img src="image_68.png" /><figcaption>Fig. 3. Current density components at $t = T/8$: (left) $J_x/J_c$; (right) $J_y/J_c$.</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>Fig. 4. Current density components at $t = T/4$: (left) $J_x/J_c$; (right) $J_y/J_c$.</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>Fig. 5. Norm $J_c^{-1}\sqrt{J_x^2 + J_y^2}$ at different instants: (left) $t = T/8$; (right) $t = T/4$.</figcaption></figure>

### C. Example with Nonuniform $n$

It is well known that in practice the $n$ exponent is not a constant: it locally varies and becomes large in regions where the temperature is close to $0K$.

In this example, we suppose that $n(r) = n_0(\exp[(4T/R)]$, with $n_0 = 1$ and $r = \sqrt{x^2 + y^2}$. We note that $n(R) = 54$ and $n(0) = 1$. We consider a transverse magnetic field with $B_0 = 15$ mT and $f = 0.5$ Hz. It is applied in the $x$-direction, and the diffusion problem of electric field is scalar.

The comparison to the results issued from *H*-formulation implemented in the Comsol software shows the validity of our approach. Fig. 6 presents the current density distribution at $t =$

<figure><img src="image_71.png" /><figcaption>Fig. 6. Current density distribution at $t = 0.5$ s, computed by DG method with $p = 4$ (left) and *H*-formulation implemented in the Comsol software (right).</figcaption></figure>

$T/4$. Near the border $n$ is large and $J/J_c$ is close to 1. When approaching the center of the cylinder, $n$ decreases and $J/J_c$ becomes smaller.

## VI. CONCLUSION

In this paper, we presented a Discontinuous Galerkin method for solving the nonlinear diffusion problems describing the evolution of the electric field in superconducting materials. The robustness of the approach was highlighted on examples where the exponent in the power constitutive law is large and uniform (for which standard finite-element codes fail to converge), and when the exponent is space-dependent.

The implementation of a semi-implicit 3D scheme based on the same approach as the one used in 2D is currently in progress.

## ACKNOWLEDGMENT

This work has been supported in part by the Belgian Science Policy (IAP P6/21) and the Belgian French Community (ARC 09/14-02).

## REFERENCES

[1] C. P. Bean, "Magnetization of hard superconductors," Phys. Rev. Lett., vol. 8, no. 6, pp. 250-253, 1962.

[2] R. Brambilla, F. Grilli, and L. Martini, "Development of an edge-element model for AC losses computation of high-temperature superconductors," Supercond. Sci. Technol., vol. 20, no. 1, pp. 16-24, 2007.

[3] A. Kameni, S. Mezani, F. Sirois, D. Netter, J. Leveque, and B. Douine, "A 2-D robust FE-FV mixed method to handle strong nonlinearities in superconductors," IEEE Trans. Magn., vol. 46, no. 8, pp. 3445-3448, 2010.

[4] F. Sirois and F. Roy, "Computation of 2-D current distribution in superconductors of arbitrary shapes using a new semi-analytical method," IEEE Trans. Appl. Supercond., vol. 17, no. 3, pp. 3836-3845, 2007.

[5] Stavrev, S. Grilli, F. Dutoit, B. Nibbio, N. Vinot, E. Klutsch, I. Meunier, G. Tixador, P. Yifeng Yang, and E. Martinez, "Comparison of numerical methods for modeling of superconductors," IEEE Trans. Magn., vol. 38, no. 2, pp. 849-852, 2002.

[6] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, "Unified analysis of discontinuous Galerkin methods for elliptic problems," SIAM J. Numer. Anal., vol. 39, no. 5, pp. 1749-1779, 2002.

[7] C. Geuzaine and J.-F. Remacle, "GMSH: A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities," Int. J. Numer. Methods Eng., vol. 79, no. 11, pp. 1309-1331, 2009.

[8] W. J. Carr, Jr., AC Loss and Macroscopic Theory of Superconductors, 2nd ed. London, U.K.: Taylor and Francis, 2001.
# A 3-D Semi-Implicit Method for Computing the Current Density in Bulk Superconductors

Abelin Kameni<fnref n="1" />, Mohamed Boubekeur<fnref n="1" />, Lotfi Alloui<fnref n="1" />, Frédéric Bouillault<fnref n="1" />, Jonathan Lambretchs<fnref n="2" />, and Christophe Geuzaine<fnref n="3" />

A semi-implicit approach is proposed for computing the current density in superconductors characterized by nonlinear vectorial power law. A nodal discontinuous Galerkin method is adopted for the spatial discretization of the nonlinear system satisfied by the components of the electric field. Explicit developments are used to construct boundary conditions to avoid the modeling of a volume around the superconducting sample. A modified Newton iterative method is introduced for solving the discrete system. Numerical examples on a 2-D superconducting plate and a 3-D superconducting cube are computed. Distributions of a component of the current density are presented and differences in the diffusive process are highlighted. The penetration time and losses are compared with those obtained with an A-V formulation solved by a finite-volume method.

**Index Terms**—Discontinuous Galerkin method, FEM, high-temperature superconductors, Maxwell's equations, Newton iterative method, semi-implicit approach.

## I. INTRODUCTION

SUPERCONDUCTING materials are increasingly introduced in electrical engineering and development of efficient numerical tools is a crucial issue for their characterization. Linear models based on critical state suggested by Bean have been widely used to evaluate losses and magnetization of low-temperature superconductors [1], but they are not suitable for high-temperature superconductors. Indeed, numerical simulations are notoriously difficult to carry out when a power law description of their behavior is assumed by J-E or E-J relations given by

$$ \frac{\vec{E}}{E_c} = \left\| \frac{\vec{J}}{J_c} \right\|^{n-1} \frac{\vec{J}}{J_c} \quad (1) $$

$$ \frac{\vec{J}}{J_c} = \left\| \frac{\vec{E}}{E_c} \right\|^{\frac{1}{n}-1} \frac{\vec{E}}{E_c} \quad (2) $$

with $\vec{J}$ the current density, $\vec{E}$ the electric field, $E_c$ the critical electric field, $J_c$ the critical current density, and the exponent $n > 1$ that characterizes the flux creep phenomena [2] ($n = 1$ corresponds to a normal conductor and $n = +\infty$ coincides with the Bean critical state).

Different discretization techniques have been implemented for computing 2-D induced fields in superconducting materials [3], [4]. For the 3-D case, the classical FEM is often used to deal with the E-J power law [6], [7]. Discretization concerns spatial operators with a nonlinear term of type

$$ \vec{\nabla} \times \vec{\nabla} (\| (\vec{J}/J_c) \|^{n-1} (\vec{J}/J_c)) \text{ or} \\ \vec{\nabla} \cdot \vec{\nabla} (\| (\vec{J}/J_c) \|^{n-1} (\vec{J}/J_c)). $$

Manuscript received June 29, 2013; revised September 13, 2013; accepted October 4, 2013. Date of current version February 21, 2014. Corresponding author: A. Kameni (e-mail: abelin.kameni@lgep.supelec.fr).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TMAG.2013.2284934

Unfortunately, during computation of the discrete solution, spurious oscillations often appear. In addition, a surrounding vacuum volume is added to set the boundary conditions. This leads to solving additional equations and increases the computational cost.

In this paper, the J-E form (2) of the constitutive law is considered and the differential system satisfied by the electric field is defined. The spatial operator is kept linear and the nonlinearity is then carried by the temporal term. This paper is based on the discretization of the scalar Laplacian operator using a nodal discontinuous Galerkin method. The interface flux terms are obtained from the symmetric interior penalty method [8]. A semi-implicit scheme that consists in solving the three nonlinear diffusion equations satisfied by the components of the electric field is proposed. The boundary conditions are defined thanks to explicit developments of Maxwell-Gauss and Maxwell-Faraday relations to avoid using a surrounding volume. They considered the demagnetizing field created by the sample. This paper is an extension of the method published previously to compute vectorial 2-D current density distributions in superconductors [9].

## II. DIFFERENTIAL SYSTEM

In a 3-D setting where the magnetic induction is $\vec{B} = (B_x, B_y, B_z)$, the electric field and the current density have three nonzero components. They satisfy a vectorial problem written as

$$ \mu_0 \frac{\partial \vec{J}}{\partial t} - \vec{\nabla} \vec{E} = - \vec{\nabla} \vec{\nabla} \cdot \vec{E}. \quad (3) $$

The main idea of this paper is to solve the diffusion equations of each component with an explicit coupling. We set $u_1 = E_x/E_c$, $u_2 = E_y/E_c$, $u_3 = E_z/E_c$, $v_1 = J_x/J_c$, $v_2 = J_y/J_c$, $v_3 = J_z/J_c$, $c = \mu_0 J_c/E_c$, and $(S_1, S_2, S_3) = -\vec{\nabla} \vec{\nabla} \cdot (\vec{E}/E_c)$ and define the following functions to consider

<footnote n="1">¹Laboratoire de Génie Electrique de Paris, Supelec, Université Paris Sud, Université Pierre et Marie Curie, Gif-sur-Yvette Cedex 91192, France</footnote>

<footnote n="2">²Centre for Systems Engineering and Applied Mechanics, Louvain-la-Neuve 1348, Belgium</footnote>

<footnote n="3">³University of Liège, Department of Electrical Engineering and Computer Science, Liège 4000, Belgium</footnote>
the 3-D behavior:

$$ \beta_1(u_1, u_2, u_3) = (u_1^2 + u_2^2 + u_3^2)^{\frac{1-n}{2n}} u_1 = v_1 \quad (4) $$

$$ \beta_2(u_1, u_2, u_3) = (u_1^2 + u_2^2 + u_3^2)^{\frac{1-n}{2n}} u_2 = v_2 \quad (5) $$

$$ \beta_3(u_1, u_2, u_3) = (u_1^2 + u_2^2 + u_3^2)^{\frac{1-n}{2n}} u_3 = v_3 \quad (6) $$

such that $\beta_1(0, u_2, u_3) = 0$, $\beta_2(u_1, 0, u_3) = 0$ and $\beta_3(u_1, u_2, 0) = 0$.

The system of three nonlinear diffusion equations is completed by fluxes boundary conditions, leading to the following system:

$$ (S) \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial \beta_1(u_1, u_2, u_3)}{\partial \beta_2(u_1, u_2, u_3)} - c^{-1} \Delta u_1 = S_1 \quad \text{in } \Omega \\[1.5ex] \displaystyle \frac{\partial \beta_2(u_1, u_2, u_3)}{\partial \beta_3(u_1, u_2, u_3)} - c^{-1} \Delta u_2 = S_2 \quad \text{in } \Omega \\[1.5ex] \displaystyle \frac{\partial \beta_3(u_1, u_2, u_3)}{\partial \beta_3(u_1, u_2, u_3)} - c^{-1} \Delta u_3 = S_3 \quad \text{in } \Omega \\[1.5ex] \begin{array}{l} \underset{\vec{v}}{\frac{\partial t}{\partial \vec{v}}} u_1 \cdot \underset{\vec{v}}{\frac{\partial t}{\partial \vec{v}}} = C_{b_1}(t) \quad \text{on } \partial \Omega \\[1.5ex] \underset{\vec{v}}{\frac{\partial t}{\partial \vec{v}}} u_2 \cdot \underset{\vec{v}}{\frac{\partial t}{\partial \vec{v}}} = C_{b_2}(t) \quad \text{on } \partial \Omega \\[1.5ex] \underset{\vec{v}}{\frac{\partial t}{\partial \vec{v}}} u_3 \cdot \underset{\vec{v}}{\frac{\partial t}{\partial \vec{v}}} = C_{b_3}(t) \quad \text{on } \partial \Omega \end{array} \end{array} \right. \quad (7) $$

where $C_{b_{1,2,3}}$ are built from Faraday's law and the assumption that $\vec{\nabla} \cdot \vec{E} = 0$ holds outside the sample.

Existence and unicity of solutions are ensured since $\beta_{1,2,3}$ are monotonous Lipschitzian functions. The proof is obtained from the study of quasi-linear elliptic and parabolic equations [10].

### III. SEMI-IMPLICIT APPROACH AND DISCRETE SYSTEM

#### A. Semi-Implicit Coupling Approach

The semi-implicit coupling approach consists in using the previous solutions $u_1^k, u_2^k$ at time $t = t_k$ for computing the solution $u_1^{k+1}$ at time $t = t_{k+1}$ (similarly for $u_2^{k+1}$ and $u_3^{k+1}$).

We note $\beta_1(u_1^{k+1}) = \beta_1(u_1^{k+1}, u_2^k, u_3^k)$, $\beta_2(u_2^{k+1}) = \beta_2(u_1^k, u_2^{k+1}, u_3^k)$ and $\beta_3(u_3^{k+1}) = \beta_3(u_1^k, u_2^k, u_3^{k+1})$. Our coupling approach leads to these expressions of nonlinearities

$$ \beta_1(u_1^{k+1}) = \left[ (u_1^{k+1})^2 + (u_2^k)^2 + (u_3^k)^2 \right]^{\frac{1-n}{2n}} u_1^{k+1} \quad (8) $$

$$ \beta_2(u_2^{k+1}) = \left[ (u_1^k)^2 + (u_2^{k+1})^2 + (u_3^k)^2 \right]^{\frac{1-n}{2n}} u_2^{k+1} \quad (9) $$

$$ \beta_3(u_3^{k+1}) = \left[ (u_1^k)^2 + (u_2^k)^2 + (u_3^{k+1})^2 \right]^{\frac{1-n}{2n}} u_3^{k+1} \quad (10) $$

#### B. Boundary Conditions

The main difficulty in superconductor modeling is to consider the demagnetizing effects. It is common to add a volume around the sample to set boundary conditions at its borders [11], but this results in increasing the computational cost. Equivalent flux boundary conditions based on computing $\vec{\nabla} u_{1,2,3}$ are proposed to avoid this technique. Their components are expressed because of explicit developments of $\vec{\nabla} \cdot \vec{E} = 0$ and $\vec{\nabla} \times \vec{E} = -\partial \vec{B} / \partial t$ on $\partial \Omega$.

These developments result in

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u_1^{k+1}}{\partial u_1^k} = - \frac{\partial u_2^k}{\partial y} - \frac{\partial u_3^k}{\partial z} \\[1.5ex] \displaystyle \frac{\partial u_1^{k+1}}{\partial y} = \frac{\partial u_2^k}{\partial x} + E_c^{-1} \frac{\partial B_z}{\partial t} \\[1.5ex] \displaystyle \frac{\partial u_1^{k+1}}{\partial z} = \frac{\partial u_3^k}{\partial x} - E_c^{-1} \frac{\partial B_y}{\partial t} \end{array} \right. \quad (11) $$

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u_2^{k+1}}{\partial u_2^k} = \frac{\partial u_1^k}{\partial y} - E_c^{-1} \frac{\partial B_z}{\partial t} \\[1.5ex] \displaystyle \frac{\partial u_2^{k+1}}{\partial y} = - \frac{\partial u_1^k}{\partial x} - \frac{\partial u_3^k}{\partial z} \\[1.5ex] \displaystyle \frac{\partial u_2^{k+1}}{\partial z} = \frac{\partial u_3^k}{\partial y} + E_c^{-1} \frac{\partial B_x}{\partial t} \end{array} \right. \quad (12) $$

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u_3^{k+1}}{\partial u_3^k} = \frac{\partial u_1^k}{\partial z} + E_c^{-1} \frac{\partial B_y}{\partial t} \\[1.5ex] \displaystyle \frac{\partial u_3^{k+1}}{\partial y} = \frac{\partial u_2^k}{\partial z} - E_c^{-1} \frac{\partial B_x}{\partial t} \\[1.5ex] \displaystyle \frac{\partial u_3^{k+1}}{\partial z} = - \frac{\partial u_1^k}{\partial x} - \frac{\partial u_2^k}{\partial y} \end{array} \right. \quad (13) $$

At $t = t_{k+1}$, the total magnetic field is defined as: $\vec{B} = \vec{B}_d + \vec{B}_d$, with $\vec{B}_d$ the uniform applied magnetic field at $t = t_{k+1}$ and $\vec{B}_d$ the demagnetized field created by the sample computed at $t_k$. The quasi-static approximation of electromagnetic fields defines the average demagnetizing field as: $\vec{B}_d = -\mu_0 \vec{\bar{N}} \vec{M}$, where $\vec{\bar{N}}$ is a tensor of the demagnetizing factors and $\vec{M}$ is the magnetization given at $t = t_k$ by: $\vec{M}^k = 0.5V \int_V \vec{r} \times \vec{J}^k dV$, with $\vec{r} = (x, y, z)$ and $V$ the volume of the sample. For the most common geometries, $\vec{\bar{N}}$ is computed from analytical formulas [12].

#### C. Discrete System

Spatial discretization of the Laplacian operator in system (7) is performed. The generated terms and the boundary conditions are represented by $F_{1,2,3}$. The nonlinear algebraic system is given by

$$ \left\{ \begin{array}{l} M_{\mathrm{DG}} \frac{\beta_1(u_1^{k+1}) - \beta_1(u_1^k)}{\beta_2(u_2^{k+1}) - \beta_2(u_2^k)} = F_1(u_1^{k+1}, u_2^k, u_3^k) \\[1.5ex] M_{\mathrm{DG}} \frac{\beta_3(u_3^{k+1}) - \beta_3(u_3^k)}{\beta_2(u_2^{k+1}) - \beta_2(u_2^k)} = F_2(u_1^k, u_2^{k+1}, u_3^k) \\[1.5ex] M_{\mathrm{DG}} \frac{\beta_3(u_3^{k+1}) - \beta_3(u_3^k)}{\partial t} = F_3(u_1^k, u_2^k, u_3^{k+1}) \end{array} \right. \quad (14) $$

where $M_{\mathrm{DG}}$ is the global mass matrix. To avoid difficulties in linearization—which introduces the fact that $\beta_1$, $\beta_2$, and $\beta_3$ are not continuous derivables at zero—a changing of discrete unknowns is proposed. The components of current density $v_1$, $v_2$, and $v_3$ are chosen as the new unknowns. They are computed from an equivalent algebraic system given by

$$ \left\{ \begin{array}{l} M_{\mathrm{DG}} \frac{v_1^{k+1} - v_1^k}{v_2^{k+1} - v_2^k} = G_1(v_1^{k+1}, v_2^k, v_3^k) \\[1.5ex] M_{\mathrm{DG}} \frac{v_2^{k+1} - v_2^k}{v_3^{k+1} - v_3^k} = G_2(v_1^k, v_2^{k+1}, v_3^k) \\[1.5ex] M_{\mathrm{DG}} \frac{v_3^{k+1} - v_3^k}{\partial t} = G_3(v_1^k, v_2^k, v_3^{k+1}) \end{array} \right. \quad (15) $$

where $G(v_1, v_2, v_3) = F(u_1, u_2, u_3)$.
<figure><img src="image_72.png" /><figcaption>Fig. 1. Partial penetration: 2-D isovalue distribution of $J_y/J_c$ at $t = 10$ ms.</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>Fig. 2. Partial penetration: 3-D isovalue distribution of $J_y/J_c$ at $t = 7.5$ ms.</figcaption></figure>

The linearization of $G$ is ensured if the inverse functions $\beta_1^{-1}$, $\beta_2^{-1}$, $\beta_3^{-1}$ exist and are derivables. Unfortunately, these derivatives are nontrivial. In accordance with the E-J power law (1), we suppose that

$$u_1^{k+1} = \beta_1^{-1}(v_1^{k+1}) = \left[ (v_1^{k+1})^2 + (v_2^k)^2 + (v_3^k)^2 \right]^{\frac{n-1}{2}} v_1^{k+1} \quad (16)$$

$$u_2^{k+1} = \beta_2^{-1}(v_2^{k+1}) = \left[ (v_1^k)^2 + (v_2^{k+1})^2 + (v_3^k)^2 \right]^{\frac{n-1}{2}} v_2^{k+1} \quad (17)$$

$$u_3^{k+1} = \beta_3^{-1}(v_3^{k+1}) = \left[ (v_1^k)^2 + (v_2^k)^2 + (v_3^{k+1})^2 \right]^{\frac{n-1}{2}} v_3^{k+1} \quad (18)$$

These inverse functions are continuous and derivable and a Newton iterative method is applied for computing $v_1^{k+1}$, $v_2^{k+1}$, and $v_3^{k+1}$ (see [9]). The time step of the numerical scheme has the same order as that is given by the CFL condition of the linear diffusion equation with $n = 1$.

## IV. NUMERICAL RESULT

We consider a superconductor sample of characteristics: $J_c = 500$ A/mm², $E_c = 10^{-7}$ V/mm, and $n = 13$. It is subjected to an external magnetic field in the z-direction: $B_a(t) = a_B t$, with $a_B = 50$ T/s. A fast ramping is considered to show robustness of this approach. The Newton iterative method needed less than five iterations at each time step $\delta t$. Three configurations are studied:

1) A 2-D superconducting plate of edge length $a = 2$ mm. The mesh has 117 000 nodes, $\delta t = 25 \mu s$. The computation time is 1 h on one core;

<figure><img src="image_74.png" /><figcaption>Fig. 3. Partial penetration: 3-D isovalue distribution of $J_y/J_c$ at $t = 11$ ms obtained with a surrounding volume.</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>Fig. 4. Full penetration: 2-D isovalue distribution of $J_y/J_c$ at $t = 19$ ms.</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>Fig. 5. Full penetration: 3-D isovalue distribution of $J_y/J_c$ at $t = 12$ ms.</figcaption></figure>

2) A 3-D superconducting cube of edge length $a = 2$ mm, whose demagnetizing factors are $N_x = N_y = N_z = 1/3$. The mesh has 80 000 nodes, $\delta t = 10 \mu s$. The computation time is 3 h on 10 cores;

3) The same 3-D superconducting cube with a surrounding cube of edge length $4a/3$ and characteristics $c = \mu_0 J_c / E_c = 0.5$, $n = 1$. The mesh has 170 000 nodes, $\delta t = 10 \mu s$. The computation time is 18 h on the same 10 cores machine. The magnetic field is imposed at border of the surrounding volume through flux expressions previously described. This results in a time shift of the penetration.

The components $J_x/J_c$ and $J_y/J_c$ have the expected symmetries due to similarities of their boundary conditions. Isovalue distributions of $J_y/J_c$ are shown in the following figures to describe partial and full penetration. They show 2-D-3-D differences in diffusive process.
<figure><img src="image_77.png" /><figcaption>Fig. 6. Full penetration: 3-D isovalue distribution of $J_y/J_c$ at $t = 13$ ms obtained with a surrounding volume.</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>Fig. 7. 3-D isovalue distribution of $J_y/J_c$ inside the surrounding conductor when the full penetration of superconducting cube is reached.</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>Fig. 8. Comparison of the dissipated energy during [0, 25 ms].</figcaption></figure>

### A. Partial Penetration

Figs. 1 and 2 show the 2-D plate at $t = 10$ ms and the 3-D cube at $t = 7.5$ ms. The not penetrated region at the center in the 2-D plate has a square form, whereas it has a cylindrical form in the 3-D case. Fig. 3 shows the partial penetration at $t = 11$ ms for the second 3-D case. We observe similarities on isovalue distributions (Figs. 2 and 3). Note that the 3-D case is similar to the 2-D one when $\vec{\nabla} \cdot \vec{E} = 0$ without the demagnetizing field.

### B. Full Penetration

The penetration times are different. The full penetration is reached at $t = 19$ ms for the 2-D plate, at $t = 12$ ms for the 3-D cube. This discrepancy is due to the demagnetizing field and $\vec{\nabla} \cdot \vec{E} \neq 0$. In the 2-D case, all the isovalues go from the corners to the center (Fig. 4), while they go from the top and bottom sides to the center in the 3-D case (Fig. 5). The penetration time of the second 3-D case is $t = 13$ ms and similarities in both 3-D distributions are observed (Figs. 5 and 6). The isovalue distribution in the surrounding volume is shown in Fig. 7.

### C. Dissipated Energy

Dissipated energy $W = \int_0^t \int_{\Omega} \vec{E} \cdot \vec{J} d\Omega dt$ is computed in different cases and compared with that obtained from a finite-volume method with A-V formulation in Fig. 8. In these published works, a surrounding volume is added and Dirichlet boundary conditions are used to impose the magnetic field without a time shift on the penetration [11]. A good agreement is observed between the results without a surrounding volume and the results obtained from an A-V formulation. When the demagnetizing effects are neglected in the 3-D case, the total magnetic field is undervalued, penetration of current density is slower, and the dissipated energy is less important in transient regime. Beyond $t = 15$ ms, these energies are close and have linear evolutions as the applied magnetic field. There is a big gap between 2-D and 3-D energies that highlights necessity to develop 3-D models.

## V. CONCLUSION

A semi-implicit approach based on a nodal discontinuous Galerkin method is proposed to describe the current density in superconducting materials. The results have been compared with those obtained from a method using a surrounding vacuum volume.

## REFERENCES

[1] A. Bossavit, “Numerical modeling of superconductors in three dimensions: A model and a finite element method,” *IEEE Trans. Magn.*, vol. 30, no. 5, pp. 3363–3366, Sep. 1994.

[2] E. Zeldov, N. M. Amer, G. Koren, A. Gupta, M. W. McElfresh, and R. J. Gambino, “Flux creep characteristics in high temperature superconductors,” *Appl. Phys. Lett.*, vol. 56, no. 7, pp. 680–683, Feb. 1990.

[3] G. Meunier and P. Tixador, “Differents formulation to models superconductors,” *IEEE Trans. Magn.*, vol. 36, no. 4, pp. 3445–3448, Jul. 2000.

[4] A. Kameni, S. Mezani, F. Sirois, D. Netter, J. Leveque, and B. Douine, “A 2-D robust FE-FV mixed method to handle strong nonlinearities in superconductors,” *IEEE Trans. Magn.*, vol. 46, no. 8, pp. 3445–3448, Aug. 2010.

[5] Z. Hong, A. M. Campbell, and T. A. Coombs, “Computer modeling of magnetisation in high temperature bulk superconductors,” *IEEE Trans. Magn.*, vol. 17, no. 2, pp. 3761–3764, Jun. 2007.

[6] F. Grilli, S. Stavrev, Y. Le Floch, M. Costa Bouzo, E. Vinot, I. Klutsch, et al., “Finite-element method modeling of superconductors: From 2-D to 3-D,” *IEEE Trans. Appl. Supercond.*, vol. 15, no. 1, pp. 17–25, Jan. 2004.

[7] M. Zhang and T. A. Coombs, “3D modeling of high-Tc superconductors by finite element software,” *Supercond. Sci. Tech.*, vol. 25, no. 1, pp. 15009–15016, 2012.

[8] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, “Unified analysis of discontinuous Galerkin methods for elliptic problems,” *SIAM J. Numer. Anal.*, vol. 39, no. 5, pp. 1749–1779, 2002.

[9] A. Kameni, J. Lambrechts, J. Remacle, S. Mezani, F. Bouillault, and C. Geuzaine, “Discontinuous Galerkin methods for computing induced fields in superconductors,” *IEEE Trans. Magn.*, vol. 48, no. 2, pp. 3445–3448, Feb. 2012.

[10] H. W. Alt and S. Luckhaus, “Quasilinear elliptic-parabolic differential equations,” *Math. Zeit.*, vol. 83, no. 6, pp. 311–341, 1983.

[11] L. Allaoui and F. Bouillault, “Numerical study of the influence of flux creep and of thermal effect on dynamic behaviour of magnetic levitation systems with a high-Tc superconductor using control volume method,” *Eur. Phys. J. Appl. Phys.*, vol. 25, no. 1, pp. 1–9, 2012.

[12] A. Aharoni, “Demagnetizing factor for rectangular ferromagnetic prisms,” *J. Appl. Phys.*, vol. 83, no. 6, pp. 3432-1–3432-3, 1998.
Formulation en champ magnétique

La formulation magnétique est celle généralement résolue par la méthode de éléments finis pour la modélisation des supraconducteurs. Avec cette formulation, on a besoin d'un volume d'air autour du matériau pour une meilleure prise en compte des effets démagnétisants. Contrairement à la méthode des éléments finis d'arêtes où le problème à résoudre a une seule inconnue vectorielle, la résolution de cette formulation par une méthode Galerkin discontinue porte sur un système de trois équations dont les inconnues sont les composantes du champ magnétique. Pour établir les équations satisfaites par les composantes ($H_1, H_2, H_3$) du champ magnétique, nous adoptons une écriture sous forme conservative de l'opérateur différentiel $\nabla \times \nabla$. On remplace alors $\nabla \times (\rho \nabla \times \mathbf{H})$ par le vecteur $(\mathrm{div}\mathbf{F}_1, \mathrm{div}\mathbf{F}_2, \mathrm{div}\mathbf{F}_3)^\mathrm{T}$ où $\mathbf{F}_1 = (0, F_z, -F_y)^\mathrm{T}$, $\mathbf{F}_2 = (-F_z, 0, -F_x)^\mathrm{T}$ et $\mathbf{F}_3 = (F_y, -F_x, 0)^\mathrm{T}$ et $(F_x, F_y, F_z)^\mathrm{T} = \rho \nabla \times \mathbf{H}$. Avec cette notation, chaque composante $H_i$ de $\mathbf{H}$ satisfait à une équation de la forme :

$$ \mu_0 \frac{\partial H_i}{\partial t} - \mathrm{div} \mathbf{F}_i = 0 \qquad (3.16) $$

La discrétisation par la méthode de Galerkin discontinue s'applique sur la formulation faible établie sur chaque élément $K$ du maillage :

$$ \mu_0 \int_K \frac{\partial H_i}{\partial t} \varphi dK - \int_K \mathbf{F}_i \cdot \nabla \varphi dK + \int_\Gamma \mathbf{F}_i \cdot \nabla \varphi d\Gamma = 0 \qquad (3.17) $$

avec $\varphi$ une fonction test et $\nu$ la normale extérieure du bord $\Gamma = \partial K$. De la même manière que dans les sections précédentes, le terme de bord est remplacé par une expression équivalente qui tient compte de la solution sur les éléments voisins. Dans ce cas, on utilise les expressions de flux construites à partir de la méthode de pénalité intérieur symétrique (Symmetric Interior Penalty method - SIP) :

$$ \int_{\Gamma} \mathbf{F}_{i} \cdot \nu_{\Gamma} \varphi d\Gamma \rightarrow \int_{\Gamma} [\mathbf{F}_{i} \cdot \nu] \{\varphi\} d\Gamma + \int_{\Gamma} \{\mathbf{F}_{i} \cdot \nu\} \cdot [\varphi] d\Gamma + I_{p_i} \qquad (3.18) $$

avec $I_{p_i}$ un terme de pénalité qui permet de satisfaire aux conditions de continuité qu'impose la physique du problème notamment celle de la continuité de la composante tangentielle. Dans notre travail nous avons montré que l'expression suivante du terme de pénalité (3.19) convient pour la modélisation des supraconducteurs et permet d'assurer la convergence de la méthode :

$$ I_{p_i} = - \int_{\Gamma} \delta[\mathbf{F}_{i}^{\Gamma}] [\varphi] d\Gamma \qquad (3.19) $$
avec $(F_1^\Gamma, F_2^\Gamma, F_3^\Gamma)^T = \mathbf{v} \times \mathbf{H} \times \mathbf{v}$ et $\delta = \alpha \rho_\Gamma h_\Gamma$ où, $\rho_\Gamma$ est la résistivité non-linéaire associée au bord $\Gamma$, $h_\Gamma$ un coefficient qui dépend de la taille de élément $K$ et $\alpha$ une constante positive.

Ce travail à fait l'objet de publications dont H-Formulation using the Discontinuous Galerkin Method for 3D modeling of superconductors publié dans *IEEE Transactions On Magnetics*.
# H-Formulation Using the Discontinuous Galerkin Method for the 3-D Modeling of Superconductors

L. Makong<sup>1,2</sup>, A. Kameni<sup>1</sup>, L. Queval<sup>1</sup>, F. Bouillault<sup>1</sup>, and P. Masson<sup>2</sup>

<sup>1</sup>Group of electrical engineering, UMR 8507 CNRS, Supelec, Université Paris Sud, Université Pierre et Marie Curie, 91192 Gif-sur-Yvette Cedex, France

<sup>2</sup>Department of Mechanical Engineering, University of Houston, Houston, TX 77204 USA

In this paper, a new approach to model efficiently high-temperature superconductors, whose electrical behavior is characterized by a non-linear power law, in 3-D is investigated. The commonly used non-linear H-formulation will be discretized based on the discontinuous Galerkin method. This numerical approach will likely provide fast simulations through parallel computations. Its application on simple examples, such as a superconducting cube subjected to an external magnetic field, will show robustness, convergence, and fast computations.

*Index Terms*—Discontinuous Galerkin (DG) method, H-formulation, high-temperature superconductors, Maxwell equations.

## I. INTRODUCTION

THE growing interest of high-temperature superconductivity applications, mostly involving machines operating in alternating current [1], [2], is synonymous with efficient design thus a precise evaluation of ac losses. Numerous difficulties associated with the highly non-linear behavior of high-temperature superconductors [3], complexities of both the geometry and the magnetic field configuration led to the ongoing development of 3-D numerical tools [4].

However, developed numerical methods, commonly using the finite-element method [5], [6], are not efficiently optimized for large scale superconductor systems such as multi-filamentary superconducting cables. Since both the memory and stiffness matrices scale with the meshed domain size, these methods are either restricted by the available memory size or the inversion of resulting large stiffness matrices. The discontinuous Galerkin (DG) method is thus implemented to make use of its natural parallel computations capabilities applied on both the mesh domain through partitioning and the resulting block diagonal matricial system [7].

Successful implementations of the method have been done to solve the formulation based on the electric field **E** [8], [9]. However, the use of demagnetization coefficients restricts this approach to model in 3-D simple geometries only. For instance, the numerical modeling of twisted multi-filamentary superconducting cables or tapes is not possible.

In this paper, the DG method will be implemented to solve the commonly used **H**-formulation coupled with the heat equation where **H** is the magnetic field and **T** the temperature. The numerical approach will be thoroughly described and applied to a superconducting cube subjected to an external sinusoidal magnetic field in both non-temperature- and temperature-

dependent cases. Comparisons of computed ac losses with the finite-element software in GetDP will be done.

## II. DIFFERENTIAL FORM

Assume a domain $\Omega$ consisting of a superconducting sub-domain $\Omega_s$ and a non-superconducting sub-domain $\Omega_r$ that are non-overlapping. Its discretization as $\Omega = \bigcup_{K \in \mathcal{T}_h} K$ will give $\mathcal{T}_h$ the mesh with tetrahedral or hexahedral elements $K$.

The combination of Maxwell equations, the magnetic linear constitutive law and the non-linear electrical power law characterizing the superconducting domain $\Omega_s$ will give the non-linear **H**-formulation. Moreover, the heat equation will be coupled to the **H**-formulation in order to have a temperature-dependent problem. The resulting differential form is

$$ egin{cases} \frac{\partial \mathbf{H}}{\partial t} + \nabla \times (\kappa(T) \nabla \times \mathbf{H}) = \mathbf{0} \quad \text{in } \Omega \\ \rho c_p(T) \frac{\partial T}{\partial t} - \nabla \cdot (\lambda(T) \nabla T) = \mathbf{J} \cdot (\sigma(T)^{-1} \mathbf{J}) \end{cases} \quad (1) $$

with $\rho$, $c_p(T)$, and $\lambda(T)$ are the mass density, the specific heat, and the thermal conductivity. The diffusivity term $\kappa(T) = 1/\sigma(T)\mu_0$ is dependent on the global domain electric conductivity $\sigma(T)$ which is either constant in the non-superconducting domain $\Omega_r$ or non-linear in the superconducting domain $\Omega_s$. The latter is expressed using the power law

$$ \mathbf{J} = \sigma(T)\mathbf{E}, \quad \text{with} \quad \sigma(T)|_{\Omega_s} = \frac{J_c(T)}{E_c} \left\| \frac{\mathbf{E}}{E_c} \right\|^{1-\frac{1}{n}} \quad (2) $$

with $J_c(T)$, $E_c$, and $n$ the critical current density, the critical electric field, and the power law index of the superconducting sub-domain $\Omega_s$. The temperature dependence of the critical current density is defined as

$$ J_c(T) = J_{c0} \frac{T_c - T}{T_c - T_0} \quad (3) $$

with $T_0$, $J_{c0}$, and $T_c$ the initial temperature of the global domain $\Omega$, the critical current density at the initial temperature, and the critical temperature of the superconducting sub-domain $\Omega_s$.

Manuscript received June 29, 2017; revised August 30, 2017; accepted October 17, 2017. Corresponding author: L. Makong (e-mail: lmakonghellnkatack@uh.edu).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TMAG.2017.2768900
Here, we applied an external magnetic field **H**<sub>a</sub> and a thermal insulation on the global domain boundary ∂Ω = Γ. Thus, we have

$$
\left\{
\begin{array}{l}
\mathbf{H} \times \mathbf{n} = \mathbf{H}_a \times \mathbf{n} \quad \text{on } \Gamma \\
\\
\nabla T \cdot \mathbf{n} = 0 \quad \text{on } \Gamma.
\end{array}
\right.
\qquad (4)
$$

III. DISCRETE VARIATIONAL FORMULATION

According to the nodal DG approach, both unknown mag-
netic field **H** and temperature **T** will be approximated as **u** and
**v** over each finite element **K**. With ℙ<sup>*m*</sup> the set of polynomials
of degree at least equal at *m*, we have

$$
\mathbf{W} = \left\{ \mathbf{w} \in (\mathbb{L}^2(\Omega))^3 : \mathbf{w}|_K \in (\mathbb{P}^m(K))^3, K \in \mathcal{T}_h \right\} \quad (5a)
$$

$$
\mathbf{V} = \left\{ v \in \mathbf{L}^2(\Omega) : v|_K \in \mathbb{F}^m(K), K \in \mathcal{T}_h \right\}. \quad (5b)
$$

Thus, the discrete variational formulation will consist in
finding on each element K, (u, v) ∈ W × V such that

$$
\mathbf{H} : \int_K \partial_t \mathbf{u} \cdot \boldsymbol{\varphi} \, dK + \int_K \kappa (\nabla \times \mathbf{u} \cdot \nabla \times \boldsymbol{\varphi}) \, dK \\
- \int_{\partial K} ((\kappa \nabla \times \mathbf{u}) \times \mathbf{n}) \cdot \boldsymbol{\varphi} \, dA = 0 \quad (6a)
$$

$$
T : \int_K \rho c_p \partial_t v \cdot \phi dK - \int_K \lambda \nabla v \cdot \nabla \phi dK \\
- \int_{\partial K} (\lambda \nabla v \cdot \mathbf{n}) \phi dA = \int_K \mathbf{j} \cdot (\sigma(v)^{-1} \mathbf{j}) \phi dK \quad (6b)
$$

with **φ** and **φ** the test functions, **j** = ∇ × **u**, and the normal **n**
is defined at each interface.

Let's assume **F**<sub>1</sub> = (0, F<sub>z</sub>, -F<sub>y</sub>)<sup>T</sup>, **F**<sub>2</sub> = (-F<sub>z</sub>, 0, F<sub>x</sub>)<sup>T</sup>,
and **F**<sub>3</sub> = (F<sub>y</sub>, -F<sub>x</sub>, 0)<sup>T</sup> with κ∇ × **u** = (F<sub>x</sub>, F<sub>y</sub>, F<sub>z</sub>)<sup>T</sup>.

The conservative form of (6a) expressed by changing the
**curl** – **curl** operator with the **div** operator will give (with
i = 1, 2, and 3) three equations of unknowns **u** = (u<sub>1</sub>, u<sub>2</sub>, u<sub>3</sub>)
such as

$$
\int_K \partial_t u_i \cdot \boldsymbol{\varphi} dK + \int_K \mathbf{F}_i \cdot \nabla \boldsymbol{\varphi} dK - \int_{\partial K} (\mathbf{F}_i \cdot \mathbf{n}) \varphi_i dA = 0 \quad (7)
$$

with

$$
\left\{
\begin{array}{l}
\nabla \times (\kappa \nabla \times \mathbf{u}) = (\nabla \cdot \mathbf{F}_1, \nabla \cdot \mathbf{F}_2, \nabla \cdot \mathbf{F}_3)^T \\
(\kappa \nabla \times \mathbf{u}) \times \mathbf{n} = (\mathbf{F}_1 \cdot \mathbf{n}, \mathbf{F}_2 \cdot \mathbf{n}, \mathbf{F}_3 \cdot \mathbf{n})^T.
\end{array}
\right.
\qquad (8)
$$

However, the continuity of the tangential components of the
magnetic field **H** has not been implemented yet. According
to the symmetric interior penalty method, equivalent interface
terms *I*<sub>*h,q*</sub> and *I*<sub>*h,q*</sub><sup>*T*</sup>, for the magnetic and thermal problem,
respectively, expressed as numerical fluxes will replace actual
interface terms in order to ensure the problem convergence.
The penalty term of *I*<sub>*h,q*</sub> will be defined such that the con-
tinuity of the tangential component of **H** is respected. The
penalty term of *I*<sub>*h,q*</sub><sup>*T*</sup> will be defined such that the continuity
of *v* and ∇*v* is respected.

IV. INTERFACE TERMS BASED ON NUMERICAL FLUXES

A. Magnetic Problem

In the magnetic problem, the equivalent interface term *I*<sub>*h,q*</sub>
will be the one proposed by Grote et al. [10] which is based
on the symmetric interior penalty method. On each face *f*
belonging to two neighboring elements *K* and *K'* or to the
global boundary *Γ*, we have the following expression:

$$
- \int_f ( [\boldsymbol{\varphi} \times \mathbf{n}] \cdot \{\kappa \nabla \times \mathbf{u}\} + [\mathbf{u} \times \mathbf{n}] \cdot \{\kappa \nabla \times \boldsymbol{\varphi}\} ) dA + I^p. \quad (9)
$$

The penalty term *I*<sup>*p*</sup> is expressed as follows:

$$
I^p = \int_f a \llbracket \mathbf{u} \times \mathbf{n} \rrbracket \cdot \llbracket \boldsymbol{\varphi} \times \mathbf{n} \rrbracket dA \quad (10)
$$

with the penalty coefficient a dependent on both the size and
diffusivity of both neighboring elements K and K'.

The quantities [[**u** × **n**]] and {**u** × **n**} represent the jump and average of the tangential components of the field **u** across each face *f* common to the neighboring elements *K* and *K'* [10]. They are defined as follows:

$$
\llbracket \mathbf{u} \times \mathbf{n} \rrbracket = \mathbf{u}_K \times \mathbf{n} - \mathbf{u}_{K'} \times \mathbf{n} \quad (11)
$$

and

$$
\{\mathbf{u} \times \mathbf{n}\} = \frac{\mathbf{u}_K \times \mathbf{n} + \mathbf{u}_{K'} \times \mathbf{n}}{2}. \quad (12)
$$

In this paper, the interface term *I*<sub>*h*</sub> is replaced by an
equivalent numerical flux *I*<sub>*h,q*</sub> given by (15), where the penalty
term *Ip* given by (10) is rewritten using the mixed product
invariance property, **a** · (**b** × **c**) = **b** · (**c** × **a**) = **c** · (**a** × **b**), in
order to respect the variational formulation based on the *div*
operator

$$
I_{h,q} = \int_f ( [\phi_i] \{\mathbf{F}_i \cdot \mathbf{n}\} + [\mathbf{F}_i \cdot \mathbf{n}] \{\phi_i\} ) dA + I^p. \quad (13)
$$

The penalty term $I^p$ becomes

$$
I^p = - \int_f a \llbracket \boldsymbol{\varphi} \rrbracket \cdot \llbracket \mathbf{n} \times \mathbf{u} \times \mathbf{n} \rrbracket dA \quad (14)
$$

with the penalized coefficient *a* > 0.

B. Thermal Problem

In the thermal problem, the equivalent interface term $I_{h,q}^T$,
based on the symmetric interior penalty method [7], will be
defined as follows:

$$
I_{h,q}^T = - \int_f ( [\lambda \nabla \phi \cdot \mathbf{n}] \{v\} - [\lambda \nabla v \cdot \mathbf{n}] \{\phi\} ) dA + I^{p,T} \quad (15)
$$

with the penalty term $I^{p,T}$ defined as

$$
I^{p,T} = \int_f b \llbracket v \rrbracket \llbracket \phi \rrbracket dA \quad (16)
$$

with penalty coefficient b > 0.
## V. BOUNDARY CONDITIONS

In order to respect the physics, Neumann and Dirichlet boundary conditions must be imposed in accordance with the numerical fluxes developed previously. The associated boundaries in both the magnetic and thermal problems, respectively, are $\Gamma^D$, $\Gamma^{D,T}$ for the Dirichlet ones and $\Gamma^N$, $\Gamma^{N,T}$ for the Neumann ones. They are defined such that $\Gamma^D \cap \Gamma^N = \emptyset$, $\Gamma^{D,T} \cap \Gamma^{N,T} = \emptyset$, $\Gamma = \Gamma^D, T \cup \Gamma^N, T$, and $\Gamma = \Gamma^D \cup \Gamma^N$. The interface common to both neighboring elements K and K' is always assumed to belong to $\Gamma^D, \Gamma^{D,T}$ or $\Gamma^N, \Gamma^{N,T}$ on the side of K.

### A. Magnetic Problem

In this problem, Dirichlet boundary condition, expressed as $\mathbf{H} = \mathbf{H}_a$ where $\mathbf{H}_a$ is the applied magnetic field, will be applied on the boundary $\Gamma^D$. Thus, the following jump term present on the penalized term of $I_{h,q}$:

$$ \llbracket \mathbf{n} \times \mathbf{u} \times \mathbf{n} \rrbracket = \frac{\mathbf{n} \times \mathbf{H}_a \times \mathbf{n} - \mathbf{n} \times \mathbf{u}_{K'} \times \mathbf{n}}{2}. \quad (17) $$

The quantities $\mathbf{u}_K$ and $\mathbf{u}_{K'}$ are the approximation of the magnetic field $\mathbf{H}$ over the neighboring elements K and K', respectively. Meanwhile, the Neumann boundary condition $\mathbf{H} \times \mathbf{n} = \mathbf{0}$, applied on $\Gamma^N$, will give the following jump term present on the penalized term:

$$ \llbracket \mathbf{n} \times \mathbf{u} \times \mathbf{n} \rrbracket = \frac{\mathbf{n} \times \mathbf{u}_{K'} \times \mathbf{n}}{2}. \quad (18) $$

### B. Thermal Problem

In this problem, Dirichlet boundary condition, expressed as $T = T_a$ where $T_a$ is the fixed temperature, imposed on the boundary $\Gamma^{D,T}$. Thus, the following jump term present on the penalized term of $I_{h,q}$:

$$ \llbracket v \rrbracket = \frac{T_a + v_{K'}^h}{2}. \quad (19) $$

The quantities $v_K$ and $v_{K'}$ are the approximation of the temperature $T$ over the neighboring elements K and K', respectively. Meanwhile, the Neumann boundary condition $\mathrm{grad}T \cdot \mathbf{n} = \mathbf{0}$, applied on $\Gamma^{N,T}$, will give the following jump term present on the penalized term:

$$ \llbracket v \rrbracket = \frac{v_{K'}}{2}. \quad (20) $$

## VI. NUMERICAL TREATMENT OF THE NON-LINEARITIES ARISING FROM THE POWER LAW

The non-linearities resulting from the superconducting sub-domain $\Omega_s$ electrical behavior of $\Omega_s$ will be dealt with explicitly.

The non-linear conductivity $\sigma^{l-1}$ of the sub-domain $\Omega_s$, evaluated at the previous time step $t_p^{l-1}$ of the problem resolution, will be used to compute the problem at following time step $t_p^l$.

<figure><img src="image_80.png" /><figcaption>Fig. 1. AC losses in the non-thermal case of a superconducting cube subjected to an external magnetic field computed with GetDP, Comsol, and the DG method with the non-linear resistivity evaluated explicitly.</figcaption></figure>

## VII. NUMERICAL RESULTS

A superconducting cube of 2 mm side length, subjected to an external magnetic field $\mathbf{H}_a = H_m \sin(2\pi ft)\mathbf{e}_y$, will be modeled using both the non-thermal coupling and thermal coupling formulation with the DG method. The mesh size was about 21 000 elements. The computed ac losses will be compared with the ones obtained with an H-formulation finite-element model implemented in GetDP and in Comsol Multi-physics. We add results obtained from Comsol Multiphysics, which is the most used software for superconductors modeling, as reference for all the validations.

### A. Non-Thermal Coupling

The superconducting behavior of the cube is characterized by a critical electric field $E_c = 10^{-7}$ V/mm, a critical current density $J_c = 100$ A/mm$^2$ and a power law exponent $n = 10$. The magnetic applied flux density amplitude is $B_m = 2$ T with a frequency $f = 50$ Hz.

Computed ac losses, Fig. 1, are equivalent. A first-order backward Euler method was used as a time integrator and a timestep of 5 µs for both approaches. More details on the computations parameters are given in Table I. Computations time are defined for a quarter of cycle. The allocated memory was about 1 GB for the new approach and about 3 GB for GetDP.

### B. Thermal Coupling

The superconducting cube used in the non-thermal coupling case will be modeled such that there is simultaneously a thermal insulation and an external magnetic field $\mathbf{H}_a = H_m \sin(2\pi ft)\mathbf{e}_y$ on the air domain boundaries. Comparisons of the computed ac losses in this thermal case will be done using the numerical approaches investigated in the non-thermal coupling case.

The superconducting behavior of the cube is characterized by a critical electric field $E_c = 10^{-7}$ V/mm, an initial critical current density $J_{c0} = 100$ A/mm$^2$, a critical temperature $T_c = 82$ K and a power law exponent $n = 10$. The magnetic
<figure><img src="image_81.png" /><figcaption>Fig. 2. Temperature distribution in Kelvins at $t = 20$ ms in the superconducting cube.</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>Fig. 3. AC losses in the thermal case of a superconducting cube subjected to an external magnetic field computed with GetDP, Comsol, and the DG method with the non-linear resistivity evaluated explicitly.</figcaption></figure>

applied flux density amplitude is $B_m = 2$ T with a frequency $f = 50$ Hz.

While in the air domain, the thermal properties are the specific heat $c_{pa} = 2040$ J/(kg.K), the mass density $\rho_a = 8.08607$ kg/m³, and the thermal conductivity $\lambda_a = 0.139$ W/(m.K), the superconducting cube is characterized by a specific heat $c_{ps} = 3$ J/(kg.K), the mass density $\rho_s = 6300$ kg/m³, and the thermal conductivity $\lambda_s = 11$ W/(m.K). The initial temperature $T_o$ is set at 77 K.

Computed ac losses, Fig. 3, remain identical. A first-order backward Euler method was used as a time integrator and a timestep of 5 µs for both approaches. More details on the computations parameters are given in Table I. The allocated memory was about 1 GB for the new approach and about 3 GB for GetDP.

The distribution, shown in Fig. 2, illustrates the elevation of temperature, below the critical state, occurring in the superconducting cube.

<table><caption>TABLE I
COMPUTATIONS PARAMETERS</caption><thead><tr><th>Case</th><th>Method</th><th>DOFs H</th><th>DOFs T</th><th>Cores</th><th>Computation time</th></tr></thead><tbody><tr><td rowspan="3">Non thermal</td><td>GetDP</td><td>23889</td><td>-</td><td>1</td><td>67 mn</td></tr><tr><td rowspan="2">DG</td><td rowspan="2">252012</td><td>-</td><td>1</td><td>128 mn</td></tr><tr><td>-</td><td>4</td><td>55 mn</td></tr><tr><td rowspan="3">Thermal</td><td>GetDP</td><td>23889</td><td>3277</td><td>1</td><td>116 mn</td></tr><tr><td rowspan="2">DG</td><td rowspan="2">252012</td><td rowspan="2">84004</td><td>1</td><td>168 mn</td></tr><tr><td>4</td><td>75 mn</td></tr><tr><td></td><td></td><td></td><td></td><td>8</td><td>51 mn</td></tr></tbody></table>

## VIII. CONCLUSION

The use of the DG method on the H-formulation successfully models high-temperature superconductors whether thermal effects are neglected or not. With parallel computations, it gives us a scalable and distributed approach in terms of both memory allocation and diagonal block matricial system resolution.

Because of the local finite-element approach used on each element of the meshed domain, both the numerical treatments of the non-linearities and the mesh domain through partitioning are distributed through parallel computations.

Unlike the DG method applied on the E-formulation [8], there is no need to compute analytically the demagnetization coefficients of the studied domain.

All those advantages make the DG method on the H-formulation suitable for the 3-D modeling of complex and large superconducting systems in real applications settings.

## REFERENCES

[1] P. J. Masson, D. S. Soban, E. Upton, J. E. Pienkos, and C. A. Luongo, “HTS motors in aircraft propulsion: Design considerations,” *IEEE Trans. Appl. Supercond.*, vol. 15, no. 2, pp. 2218–2221, Jun. 2005.

[2] S. D. Chen et al., “Design of a HTS magnet for application to resonant X-ray scattering,” *IEEE Trans. Appl. Supercond.*, vol. 21, no. 3, pp. 1661–1664, Jun. 2011.

[3] E. Zeldov, N. M. Amer, G. Koren, A. Gupta, M. W. McElfresh, and R. J. Gambino, “Flux creep characteristics in high-temperature superconductors,” *Appl. Phy. Let.*, vol. 56, no. 7, pp. 680–683, 1990.

[4] A. Bossavit, “Numerical modelling of superconductors in three dimensions: A model and a finite element method,” *IEEE Trans. Magn.*, vol. 30, no. 5, pp. 3363–3366, Sep. 1994.

[5] F. Grilli et al., “Finite-element method modeling of superconductors: From 2-D to 3-D,” *IEEE Trans. Appl. Supercond.*, vol. 15, no. 1, pp. 17–25, Mar. 2005.

[6] M. Zhang and T. A. Coombs, “3D modeling of high-Tc superconductors by finite element software,” *Supercond. Sci. Tech.*, vol. 25, no. 1, pp. 15009–15016, 2012.

[7] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, “Unified analysis of discontinuous Galerkin methods for elliptic problems,” *SIAM J. Numer. Anal.*, vol. 39, no. 5, pp. 1749–1779, 2002.

[8] A. Kameni, M. Boubekeur, L. Alloui, F. Bouillault, J. Lambrechts, C. Geuzaine, “A 3-D semi-implicit method for computing the current density in bulk superconductors,” *IEEE Trans. Magn.*, vol. 50, no. 2, Feb. 2014, Art. no. 7009204.

[9] A. Kameni, J. Lambrechts, J.-F. Remacle, S. Mezani, F. Bouillault, and C. Geuzaine, “Discontinuous Galerkin method for computing induced fields in superconducting materials,” *IEEE Trans. Magn.*, vol. 48, no. 2, pp. 591–594, Feb. 2012.

[10] M. J. Grote, A. Schneebeli, and D. Schötzau, “Interior penalty discontinuous Galerkin method for Maxwell’s equations: Energy norm error estimates,” *J. Comput. Appl. Math.*, vol. 204, no. 2, pp. 375–386, 2007.
3.4 Modélisation des câbles supraconducteurs torsadés

Les câbles supraconducteurs ont été développés pour diverses applications telles que les moteurs, les transformateurs, les limiteurs de courant, les aimants ou encore les câbles de transmission d'énergie électrique. Ceux-ci sont notamment destinés à remplacer les câbles en cuivre. Avec ces nouveaux câbles il est possible de transporter une puissance cinq fois plus grande qu'avec les câbles classiques en cuivre. Ces câbles présentent aussi l'avantage d'occuper moins d'espace et permettent de réduire considérablement les pertes d'énergie qui sont de l'ordre de 10% avec les câbles conventionnels en cuivre [86]. Ils sont généralement constitués de plusieurs dizaines de filaments supraconducteurs fins intégrés à une matrice résistive qui sert de protection aux filaments contre d'éventuelles sollicitations mécaniques. Ceux-ci sont torsadés dans le but de concentrer les courants induits dans les filaments supraconducteurs et ainsi réduire les pertes Joule qui se développent dans la matrice resistive [25]. L'évaluation précise des pertes générées par ces câbles est important pour améliorer le dimensionnement des systèmes. Les outils numériques permettant d'évaluer ces pertes doivent être robustes pour faire face à la complexité introduite par le torsadage des filaments et le comportement électrique non-linéaire du matériau.

La modélisation des câbles supraconducteurs s'est appuyée pendant plusieurs décennies sur des formules analytiques basées sur le modèle de Bean [23]. Depuis l'introduction des supraconducteurs à haute température critique dans la fabrication de ces câbles, ces travaux de modélisations s'appuient sur des méthodes numériques notamment la méthode des éléments finis [2]. La plupart des modèles numériques proposés sont bi-dimensionnels et permettent d'évaluer approximativement les grandeurs telles que les pertes ou l'aimantation de ces câbles avec une précision moindre [95, 94]. La précision de ces résultats peut être améliorée notamment avec l'utilisation de modèles tri-dimensionnels. Malheureusement, l'utilisation des méthodes de discrétisation telle que la méthode des éléments finis pour la modélisation tri-dimensionnelle de ces câbles est très coûteuse car les systèmes à résoudre sont grands, ce qui nécessite une importante puissance de calcul pour réaliser les opérations d'inversions matricielles. Dans des travaux récent, de telles simulations ont été réalisées sur des câbles supraconducteurs torsadés [37]. Ces travaux ont mis en évidence les coûts de calculs importants que requiert ce type de modélisation. Afin de réduire ces coûts de calcul, il devient indispensable de mettre en place des techniques permettant de réduire la taille des systèmes à résoudre. Nous avons réalisé des développements sur le thème de la réduction de modèle pour la modélisation de ces câbles supraconducteurs.
3.4.1 Modèle équivalent des filaments torsadés

L'idée exploitée pour construire ce modèle équivalent a été de transformer le problème initial, qui est celui d'un câble aux filaments torsadés soumis à champ magnétique de direction connue, en un problème équivalent, qui est celui d'un câble aux filaments droits soumis à un champ magnétique dont la direction subit une transformation correspondant au torsadage filaments. Cette approche permet de simplifier la construction de la géométrie, de réduire la taille du maillage et le temps de calcul. De plus, l'utilisation des filaments droits améliore le conditionnement des matrices ce qui a un impact positif sur la convergence de la résolution numérique.

Ce travail a fait l'objet de l'article présenté ci-après **3-D Equivalent Model to Compute the Electro-Magnetic Behavior of Twisted Multi-Filamentary Superconductors Wires** publié dans la revue *IEEE Transactions On Magnetics*.
# 3-D Equivalent Model to Compute the Electro-Magnetic Behavior of Twisted Multi-Filamentary Superconductors Wires

L. Makong¹,², A. Kameni¹, F. Bouillault¹, and P. Masson²

¹Group of electrical engineering-Paris, UMR 8507 CNRS, Supelec, Université Paris Sud, Université Pierre et Marie Curie, 91192 Gif-Sur-Yvette, France

²Department of Mechanical Engineering, University of Houston, Houston, TX 77204 USA

In this paper, we investigate new methods to model in 3-D, efficiently and simply, twisted multi-filamentary superconductors subjected to an external magnetic field. Several simulation cases, based on the H-formulation implemented in GetDP, will help derive through appropriate methods, such as changing the reference from where the problem is seen and described to a reference where everything becomes simple to deal with, a simple equivalent 3-D problem to solve efficiently and quickly.

**Index Terms**—Finite-element methods, high-temperature superconductors, Maxwell equations, multi-filamentary.

## I. INTRODUCTION

SEVERAL high-temperature superconductivity applications, such as motors design for aircraft propulsion [1] and magnets design for medical imaging [2], use twisted multi-filamentary superconductors wire [3] for their windings.

Both the fine superconducting filaments structure and the twisting process greatly help reduce the ac losses by keeping the filaments as much uncoupled as possible. However, a numerical and precise evaluation of ac losses generated is complex even with the use of robust 3-D numerical tools [7]–[9].

Complexities encompass the highly non-linear electrical behavior of superconductors [5], [6] and the twisted superconducting filaments geometry. In order to compute efficiently the wire behavior, its meshed domain must require a mesh refinement on the curved parts of each twisted superconducting filaments. It will add more degrees of freedom to the already large matricial system derived. In addition, the wire model is not prone to the use of symmetry boundary conditions because of the twist transformation. The latter gives dissimilar meshed surfaces at the wire’s top and bottom. The symmetry boundary conditions are essential to reduce the modeling of the complete wire. They can allow to model only a truncated part of the wire like a single twist pitch.

Thus, it is important to develop methods that simplify the problem of modeling twisted multi-filamentary superconducting wire. It will allow us to perform more in-depth design and optimization studies for real world applications such as superconducting machine design.

After a description of the problem, based on the H-formulation [8] using the finite-element method, we investigate several simulation cases. The goal is to derive an equivalent simple 3-D problem to model quickly and precisely

the electro-magnetic behavior of twisted multi-filamentary superconducting wire subjected to an external magnetic field. All the models involved in this paper are implemented in the finite-element open source code GetDP [10].

## II. PROBLEM DEFINITION AND FORMULATION

The studied domain $\Omega$ consists of a superconducting sub-domain $\Omega_s$, which are all the filaments, embedded in a non-superconducting sub-domain $\Omega_{ns}$. The non-superconducting sub-domain includes both the resistive matrix and the air domain around the entire wire.

In order to compute the electro-magnetic behavior of the studied wire, the H-formulation will be implemented numerically using the finite-element method. The resulting variational formulation will consist in finding $\mathbf{H} \in H(\text{curl}; \Omega)$ such that

$$ \int_{\Omega} \frac{\partial \mathbf{H}}{\partial t} \cdot \varphi d\Omega + \int_{\Omega} \kappa \cdot \text{curl} \mathbf{H} \cdot \text{curl} \varphi d\Omega + I^{B,D} = 0, \quad \forall \varphi \in H(\text{curl}; \Omega) \quad (1) $$

with the boundary term $I^D = \int_{\Gamma^D} \lambda \cdot \varphi dA + \int_{\Gamma^D} (\mathbf{H} - \mathbf{H}_a) \cdot \lambda' dA$. The external magnetic field $\mathbf{H}_a$ is imposed on $\Gamma^D \in \partial\Omega$ as a Dirichlet boundary condition using Lagrange coefficients $\lambda$ and $\lambda'$.

Both the unknown magnetic field $\mathbf{H}$ and the basis vector function $\varphi$ are projected along edges elements. It will ensure the continuity of the tangential component of the magnetic field.

Moreover, the non-linear electrical behavior of superconductors is characterized by a power law [6] with a non-linear resistivity in $\Omega_s$

$$ \rho = \frac{E_c}{J_c} \left\| \frac{\mathbf{J}}{J_c} \right\|^{n-1} . \quad (2) $$

The quantities $E_c$, $J_c$, and $n$ are the critical electric field, the critical current density, and the power law index. However, the resistivity $\rho$ will remain constant in both the resistive matrix and the air domain. Numerically, the associated power

Manuscript received June 27, 2017; accepted October 30, 2017. Corresponding author: L. Makong (e-mail: lmakong@hllnkatack@uh.edu).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TMAG.2017.2771438
<figure><img src="image_83.png" /><figcaption>Fig. 1. 3-D current density distribution inside the filament at $t = T/4$ s.</figcaption></figure>

law will be linearized in the variational formulation. It will allow stability and convergence of the problem using the finite-element method.

### III. EFFECTS OF A TRANSVERSE MAGNETIC FIELD ON A TWISTED MONO-FILAMENT SUPERCONDUCTING WIRE

A twist pitch of a mono-filament superconducting wire subjected to a sinusoidal transverse magnetic field is investigated. The wire, completely surrounded by an air domain, consists of a single twisted superconducting filament embedded in a Niobium matrix.

Finite-element methods' computations of this problem, using GetDP, gave an unexpected current density distribution inside the wire: three current loops locally created, as shown in Fig. 1.

This result, although not intuitively predictable, differs from the expected current distribution of a similar magnetic configuration applied to a straight mono-filament superconducting wire. One current loop will be created inside the straight filament.

Differences between a twisted and a straight filament subjected to a similar magnetic field assume that the twisted superconducting filament observes, in the Frenet frame (T, N, and K) moving along the resulting helical trajectory, a transverse magnetic field with a directional change. The resulting helical trajectory is associated with the geometric transformation of straight filaments into twisted filaments. T, N, and K are, respectively, the unit vector tangent, the normal unit vector, and the binormal unit vector.

In order to confirm the previous assumption, a mathematical expression of the transformation is derived for numerical modeling purposes. The transformation assumes a fixed Frenet frame (T, N, and K), where the twisted filament becomes a straight one.

While the initial frame (x, y, and z) is assumed fixed before the transformation, it moves, after the transformation, along the resulting helical trajectory in the Frenet frame (T, N, and K).

Such a transformation also applies to the physics of the domain. Because the new frame remains fixed, differential operators describing the domain's physics will remain unchanged. However, the external boundary conditions they are subjected to might likely change. For instance, the external magnetic field applied to the superconducting wire must be transformed in the new frame, as shown in Fig. 2.

<figure><img src="image_84.png" /><figcaption>Fig. 2. 3-D filament subjected to a transverse magnetic field in the Frenet frame.</figcaption></figure>

<table><caption>TABLE I
MODELS PARAMETERS</caption><thead><tr><th>Quantity</th><th>Value</th></tr></thead><tbody><tr><td>Filament radius (μm)</td><td>60</td></tr><tr><td>Matrix radius (mm)</td><td>0.4</td></tr><tr><td>Critical current density (A/mm²)</td><td>5</td></tr><tr><td>Power law index n</td><td>10</td></tr><tr><td>Matrix conductivity (S/m)</td><td>6·10<sup>9</sup></td></tr><tr><td>Magnetic field amplitude B<sub>max</sub> (mT)</td><td>100 or 50</td></tr><tr><td>Frequency f (Hz)</td><td>50</td></tr><tr><td>Magnetic field norm</td><td>B<sub>m</sub>sin(2πf)</td></tr></tbody></table>

Thus the mathematical expression of the geometric transformation applied to the magnetic field **H** will be the following:

$$ \mathbf{H}_{[T,N,K]} = \begin{bmatrix} -a \cdot \sin(\theta) & a \cdot \cos(\theta) & b \\ \cos(\theta) & \sin(\theta) & 0 \\ -b \cdot \sin(\theta) & b \cdot \cos(\theta) & -a \end{bmatrix} \cdot \mathbf{H}_{[x,y,z]} \quad (3) $$

with $a = r/L$, $b = c/L$, $\theta = z/c$, $c = p/2\pi$, and $L^2 = r^2 + c^2$. The helical trajectory parameters *p* and *r* are the twist pitch and the radius.

In the Frenet frame, we computed the current density distribution of the previously mentioned twisted mono-filament superconducting wire subjected to an external magnetic field. However, the applied magnetic field is now the transformed form of the external magnetic field used in the initial frame (x, y, and z) shown in the previous case. We, therefore, have an external magnetic field with a directional change along the helical trajectory associated with the transformation, as shown in Fig. 2. A single current loop appeared inside the filament. The studied case computed in the Frenet frame is thus equivalent to the case of a straight mono-filament superconducting wire subjected to a unidirectional transverse magnetic field.

### IV. VALIDATION OF THE GEOMETRIC TRANSFORMATION ON A TWISTED BI-FILAMENTS SUPERCONDUCTING WIRE SUBJECTED TO A TRANSVERSE MAGNETIC FIELD

The geometric transformation described above will be investigated for a twisted bi-filaments superconducting wire subjected to a transverse magnetic field. Both superconducting filaments are embedded in a Niobium matrix with an air domain surrounding both filaments and matrix. The model parameters are shown in Table I.

The validation goal is to compare modeling cases with and without the geometric transformation in order to check
<figure><img src="image_85.png" /></figure>

<figure><img src="image_86.png" /></figure>

<figure><img src="image_87.png" /><figcaption>Fig. 3. AC losses comparison, with $B_{max}$ = 100 mT, of both twisted and its equivalent straight bi-filamentary superconductors in the initial (in blue) and Frenet frames (in red), respectively. (a) $p$ = 2 mm. (b) $p$ = 4 mm. (c) $p$ = 10 mm.</figcaption></figure>

<figure><img src="image_88.png" /></figure>

<figure><img src="image_89.png" /></figure>

<figure><img src="image_90.png" /><figcaption>Fig. 4. AC losses comparison, with $B_{max}$ = 50 mT, of both twisted and its equivalent straight bi-filamentary superconductors in the initial (in blue) and Frenet frames (in red), respectively. (a) $p$ = 2 mm. (b) $p$ = 4 mm. (c) $p$ = 10 mm.</figcaption></figure>

any equivalence results-wise. In both cases, it will involve a parametric study with the twist pitch $p$ and the magnetic field amplitude $B_{max}$ as variables. While the twist pitch value will vary as 2, 4, and 10 mm, the magnetic field amplitude will either be 50 or 100 mT.

Comparisons of modeling cases will focus on the computed ac losses generated by the superconducting filaments. For each defined parameters $p$ and $B_{max}$, computed ac losses, in both modeling case in the initial frame and its associated transformed case in the Frenet frame, will be compared.
<figure><img src="image_91.png" /><figcaption>Fig. 5. AC losses comparison, with $B_{max}$ = 100 mT, of both twisted and its equivalent straight six-filament superconductors in the initial (in blue) and Frenet frames (in red), respectively.</figcaption></figure>

### A. Cases With $B_{max}$ = 100 mT

The superconducting filaments are fully penetrated by the current in this case. Moreover, we notice some discrepancies between ac losses in both the initial (in blue) and Frenet (in red) frames, as shown in Fig. 3. These discrepancies, though small, are mainly notable in cases where the twist pitch is small ($p = 2$ mm and $p = 4$ mm).

However, it is observed that these discrepancies are reducing with an increase of the twist pitch. Thus, we have an ac losses equivalence, for the given filament radius, in both the initial and Frenet frames when the twist pitch is 10 mm. This result is consistent with the fact that with a longer wire, the helical trajectory, described by the transformed magnetic field, is getting closer to a straight trajectory. It will thus lead to a transformed magnetic field equivalent to the transverse magnetic field in the initial frame.

### B. Cases With $B_{max}$ = 50 mT

In this paper, there is no notable discrepancies between ac losses in both the initial (in blue) and Frenet (in red) frames, as shown in Fig. 4.

We have an ac losses equivalence in both the initial and Frenet frames regardless of the twist pitch length. Moreover, the superconducting filaments are not fully penetrated.

## V. VALIDATION OF THE GEOMETRIC TRANSFORMATION ON A TWISTED MULTI-FILAMENTARY (SIX FILAMENTS) SUPERCONDUCTING WIRE SUBJECTED TO A TRANSVERSE MAGNETIC FIELD

In this case, we have a wire that includes six twisted filaments of radius $r_f = 50\mu$m embedded in a Niobium matrix of radius $r_m = 0.6$ mm. The six filaments are distributed on two layers of radius $r_{l1} = 0.1$ mm and $r_{l2} = 0.35$ mm, respectively. Two filaments belong to the interior layer while the four left are located on the exterior one.

The magnetic properties of the problem are identical with the case of the twisted bi-filaments wire. Moreover, the magnetic field amplitude is set at 100 mT with a twist pitch of 4 mm.

As shown in Fig. 5, the computed ac losses in both the initial and the Frenet frame are equivalent for the studied case.

Moreover, it only took one week to compute the problem in the Frenet frame, whereas the computations, in the initial frame, lasted two weeks.

## VI. CONCLUSION

In this paper, we investigated a geometric transformation, which simplifies the 3-D modeling of twisted multi-filamentary superconducting wires subjected to a transverse magnetic field. The transformation allows the study of the wires in the Frenet frame moving along the helical trajectory resulting from twisting filaments.

This change of frame has been extended to the case of twisted bi-filamentary superconducting. A parametric study based on the magnetic field amplitude and the twist pitch was done. We noticed an equivalence, regardless of the twist pitch, between the transformed problem and the initial one for small magnetic field. However, small discrepancies are present, especially for small twist pitch, when the magnetic field is high. In the case involving six filaments, we found a notable equivalence for a small twist pitch in high field.

This equivalence, through the described transformation, works for non-magnetic matrix only. This methodology also works for multi-filamentary wires with normal conductors as filaments. As shown in the case of a wire with six filaments, the transformed problem reduces greatly the computation time.

In this Frenet frame, the modeling of realistic twisted multi-filamentary superconducting wires will naturally be prone to the modeling of any truncated wires parts. In fact, the transformation conserves the physics along the wire due to the moving trajectory of the external magnetic field.

For the future works, the presented parametric study will be extended in order to characterize the influence of the filament radius.

## REFERENCES

[1] P. J. Masson, D. S. Soban, E. Upton, J. E. Pienkos, and C. A. Luongo, “HTS motors in aircraft propulsion: Design considerations,” *IEEE Trans. Appl. Supercond.*, vol. 15, no. 2, pp. 2218–2221, Jun. 2005.

[2] S. D. Chen et al., “Design of a HTS magnet for application to resonant X-ray scattering,” *IEEE Trans. Appl. Supercond.*, vol. 21, no. 3, pp. 1661–1664, Jun. 2013.

[3] C. J. Christopherson and G. N. Riley, “Development of twisted high-temperature superconductor composite conductors,” *Appl. Phys. Lett.*, vol. 66, no. 17, pp. 2277–2279, 1995.

[4] G. Meunier, G. Meunier, and P. Tixador, “Different formulations to model superconductors,” *IEEE Trans. Magn.*, vol. 36, no. 4, pp. 3445–3448, Jul. 2000.

[5] C. P. Bean, “Magnetization of high-field superconductors,” *Rev. Mod. Phys.*, vol. 36, no. 1, pp. 31–39, 1964.

[6] E. Zeldov, N. M. Amer, G. Koren, A. Gupta, M. W. McElfresh, and R. J. Gambino, “Flux creep characteristics in high-temperature superconductors,” *Appl. Phy. Lett.*, vol. 56, no. 7, pp. 680–683, 1990.

[7] A. Bossavit, “Numerical modelling of superconductors in three dimensions: A model and a finite element method,” *IEEE Trans. Magn.*, vol. 30, no. 5, pp. 3363–3366, Sep. 1994.

[8] F. Grilli et al., “Finite-element method modeling of superconductors: From 2-D to 3-D,” *IEEE Trans. Appl. Supercond.*, vol. 15, no. 1, pp. 17–25, Mar. 2005.

[9] M. Costa et al., “3D modeling of coupling between superconducting filaments via resistive matrix in AC magnetic field,” *IEEE Trans. Appl. Supercond.*, vol. 13, no. 2, pp. 3634–3637, Jun. 2003.

[10] P. Dular and C. Geuzaine. *GetDP: A General Environment for the Treatment of Discrete Problems*. Accessed: May 5, 2017. [Online]. Available: http://www.geuz.org/getdp
### 3.4.2 Réduction du modèle équivalent

Le modèle réduit pour le calcul de pertes dans les câbles supraconducteurs multi-filaments torsadés découle du modèle équivalent avec des filaments droits proposé au paragraphe précédent. En effet pour un câble aux filaments droits, la section est identique en tout point de la hauteur, ce qui est simplifié l'utilisation des conditions périodiques. Dans nos développements, on a établi que, pour un câble dont la longueur est suffisamment grande, il suffit de déterminer les pertes sur une fine tranche de ce câble pour connaitre les pertes sur l'ensemble du câble. Le modèle réduit de calcul de pertes dans les câbles supraconducteurs multi-filamentaires torsadés est défini sur une fine tranche de câble aux filaments droits, sur laquelle on impose des conditions de périodicité à ses extrémités. Ce travail et sa validation par comparaison à des mesures expérimentales ont été récemment proposés dans l'article ci dessous **Reduced model to computed losses of Twisted Multi-Filamentary Superconductors Wires** soumis et accepté pour publication dans la revue *IEEE Applied Superconductivity*.
# Reduced model to compute AC losses of twisted multi-filamentary superconductors

A. Kameni<fnref n="1" />, L. Makong<fnref n="1" />, F. Bouillault<fnref n="1" />, P. J. Masson<fnref n="2" />

e-mails: abelin.kameni@geeps.centralesupelec.fr

**Abstract**—Classical numerical techniques remain very computationally expensive to model the magnetization of twisted multi-filamentary superconductors in 3D because of the required large system of equations to solve. In this paper a reduced model is proposed to approximate the AC losses of these types of wire. It is established through the projection of the physical problem in a helical Frenet basis and the use of periodic boundary conditions on a slice of the wire. The reduced model on the slice leads to systems of more reasonable size leading to smaller and faster simulations than when the whole wire is modeled. Simulation results show that the proposed method leads to computed losses similar to those obtained on a whole wire. Several simulation cases are proposed to prove the validity of the proposed model, and a comparison to some AC losses measurements is proposed to show the effectiveness of the proposed approach. The implementation relies on the H-formulation and is solved with the open-source solver Getdp.

**Index Terms**—Finite element method, H-formulation, modeling, AC losses, Twisted multi-filamentary superconductors.

## I. INTRODUCTION

DURING the last two decades, superconducting materials are increasingly introduced in large scale applications such as, for example, magnetic resonance imaging, fault current limiter, superconducting motors and generators. Superconducting tapes and wires are often used to generate high magnetic fields without the need for magnetic cores in many of these applications, thus leading to more compact and efficient systems [1]. In the case of round wires, they usually consist of fine superconducting filaments embedded in a resistive matrix to protect them against mechanical strain and provide electro-thermal stability. The configuration of the wire is so that the transport current is carried entirely by the superconducting filaments while the resistive matrix only sees transient currents due to eddy currents, coupling losses or quench. To minimize the current in the matrix material, the superconducting filaments are twisted during manufacturing. It allows for a significant reduction of the coupling losses by keeping all the filaments as much uncoupled as possible[2]. A numerical and precise evaluation of AC losses generated in such wires is important for an effective design of superconducting machines and their appropriate cryogenic systems. Unfortunately, because of the complicated geometry (aspect ratio, multiple fine filaments, etc.), the fine filaments require a very fine mesh which, coupled to the non-linearity

introduced by the power law relation that characterizes the electric behavior of the superconducting materials, reders the use of classical techniques such as finite element method challenging and very expensive in terms of resources with computation times often of several weeks. However, some work proposes the use of periodic boundary conditions to reduce the computation domain [3]. In the proposed model, the length of the wire is reduced to one sixth of the twist pitch and periodic boundary conditions are used on the geometrically congruent end faces. Unfortunately, despite this length reduction, computational times still remain significant, reaching several days per simulation.

The aim of this paper is to propose a simplification of the 3D modeling of twisted multi-filamentary superconducting wires for fast and accurate AC losses computation. In our previous published work, we have proposed an equivalent model which consists in replacing the twisted wire by a straight one [4]. The geometry of this equivalent model is easier to define and the required mesh to accurately compute the losses are coarser than for a twisted geometry. Despite such simplifications, modeling of long twisted superconducting wires remains difficult and very computationally expensive despite the use of an equivalent model based on straight filaments. This is due to the great length of wires that leads to very large matrix systems to solve.

In this paper, an approach aimed at reducing the equivalent 3D problem size is proposed. This reduced model consists in modeling a slice of the straight wire with periodic boundary conditions at its ends. This work shows that it is possible to approximate the AC losses per unit length of the whole twisted superconducting wire by computations performed on a slice of the equivalent straight wire. As the reduced model is based on straight wire, the condition of geometrically congruent end faces needed to impose the periodic boundary conditions in [3] is no longer necessary, and the length of the slice can be very small. A validation of the results obtained by comparison to measurements is proposed to show relevance and effectiveness of this work. Note that in this paper AC losses refers to magnetization losses only. The numerical results proposed in this paper are obtained with a computer Intel Core i5, 2.3 Ghz, 16 Go RAM.

<footnote n="1">¹ GeePs - Group of electrical engineering - Paris, UMR CNRS 8507, CentraleSupélec, Université Paris-Sud, Université Paris-Saclay, Sorbonne Université, 11 rue Joliot-Curie, Plateau de Moulon 91192 Gif-sur-Yvette CEDEX, France.</footnote>

<footnote n="2">² AML Superconductivity and Magnetics, 1720 Main St, Bldg 4, Palm Bay, FL-32905, USA</footnote>
## II. EQUIVALENT 3D MODEL OF TWISTED MULTI-FILAMENTARY WIRE

In previously published work [4], it has been shown that the model defined by a twisted multi-filamentary wire subjected to an external axial or transverse magnetic field (Figure.1) can be replaced by an equivalent one defined as a straight multi-filamentary wire subjected to an external magnetic field which initial direction undergoes a helical transformation such as the twist process (Figure 2). In other words, a twisted multi-filamentary wire subjected to an H[x,y,z] magnetic field in the (x, y, z) frame become a straight multi-filamentary wire subjected to an H[T,N,K] in the Frenet frame (T, N, K). In accordance to the geometrical relation between the Frenet and cartesian frames (Figure 3), the external magnetic field H[T,N,K] is obtained from the following relation:

$$ H_{[T,N,K]} = \begin{bmatrix} -a \sin \theta & a \cos \theta & b \\ \cos \theta & \sin \theta & 0 \\ -b \sin \theta & b \cos \theta & -a \end{bmatrix} \cdot H_{[x,y,z]} \quad (1) $$

with, $a = r/L$, $b = c/L$, $\theta = z/c$, $c = p/2\pi$ and $L^2 = r^2 + c^2$. The wire length is L and the helical trajectory parameters p, r are respectively the twist pitch and the radius.

<figure><img src="image_92.png" /><figcaption>Fig. 1. A twisted superconducting filament subjected to a transverse magnetic field in the initial (x, y, z) frame</figcaption></figure>

Note that, because the Frenet frame remains fixed, differential operators describing the domain's physics will remain unchanged, thus the same H-formulation is solved in both frames. The weak form of this formulation is given by:

$$ \int_{\Omega} \frac{\partial \mathbf{H}}{\partial t} \cdot \psi d\Omega + \int_{\Omega} \frac{\rho}{\mu_0} (\nabla \times \mathbf{H}) \cdot (\nabla \times \psi) d\Omega + I^{Bnd} = 0 \quad (2) $$

where, $\Omega$ is the calculation domain, $\rho$ is the resistivity, $\mu_0$ is the permeability of the vacuum, $I^{Bnd}$ is the boundary term and $\psi$ is a test function taken in the functional space $H(\text{curl}, \Omega)$. The boundary term $I^{Bnd}$ allows to impose the external magnetic field $H_a$ on the boundary $\Gamma = \partial\Omega$. It is

<figure><img src="image_93.png" /><figcaption>Fig. 2. The equivalent model of the twisted superconducting filament in the Frenet frame (T, N, K): A straight superconducting filament subjected to a helical transverse magnetic field</figcaption></figure>

<figure><img src="image_94.png" /><figcaption>Fig. 3. Geometrical relation between the Frenet frame (T, N, K) and the cartesian frame (x, y, z)</figcaption></figure>

a weak form of a Dirichlet boundary condition that uses the Lagrange coefficients $\lambda$ and $\lambda'$ defined in $H(\text{curl}, \Omega)$.

$$ I^{Bnd} = \int_{\Gamma} \lambda \cdot \psi d\Gamma + \int_{\Gamma} (\mathbf{H} - \mathbf{H}_a) \lambda' d\Gamma \quad (3) $$

The resistivity $\rho$ has a constant value except in the superconducting filaments where it is described by a nonlinear relationship:

$$ \rho = \frac{E_c}{J_c} \left\| \frac{\mathbf{J}}{J_c} \right\|^{
\begin{array}{c}
n-1 \\
n \\
n
\end{array}
} \quad (4) $$

with, $E_c = 10^{-4}$ V/m, $J_c$ the critical current density and n a power law exponent that characterizes the flux creep