Helicoidal Transformation Method for Finite Element Models of Twisted Superconductors

This article deals with the modeling of superconducting and resistive wires with a helicoidal symmetry, subjected to an external field and a transport current. Helicoidal structures are 3-D, and therefore yield computationally intensive simulations in a Cartesian coordinate system. We show in this article that by working instead with a helicoidal system of coordinates, the problem to solve can be made 2-D, drastically reducing the computational cost. We first introduce the state-of-the-art approach and apply it on the <inline-formula><tex-math notation="LaTeX">$h$</tex-math></inline-formula>-<inline-formula><tex-math notation="LaTeX">$\phi$</tex-math></inline-formula>-formulation with helicoidally symmetric boundary conditions (e.g., axial external magnetic field, with or without transport current), with an emphasis on the function space discretization. Then, we extend the approach to general boundary conditions (e.g., transverse external magnetic field), and we present numerical results with linear materials. In particular, we discuss the frequency-dependent losses in composite wires made of superconducting filaments embedded in a resistive matrix. Finally, we provide outlook to the application of the generalized model with nonlinear materials.


I. INTRODUCTION
L OW-TEMPERATURE superconducting composite wires usually consist of a large number of superconducting filaments embedded in a conducting matrix.This matrix helps in redistributing current between filaments, but has the side effect of coupling the filaments in the presence of an external transverse time-varying magnetic field.This coupling can however be reduced by twisting the composite wire [1], [2].The resulting geometry is not invariant along the wire axis and leads to a computationally intensive three-dimensional (3D) modelling [3], [4], [5], [6].
Approximate models exploiting the multifilamentary structure of this kind of wires have been investigated to reduce the computational cost, such as in [7] and [8], where coupling currents in the conducting matrix are accounted for in a 2D finite element model by introducing equivalent resistances between the filaments.Alternatively, a Frenet frame is used in [9] to simplify the definition of the 3D geometry, and AC losses are approximated by considering a fraction of the pitch length of the wire in a 3D model, or a cross section of the wire in a 2D model.Homogenization techniques involving anisotropic materials have also been considered [3].Finally, parallelization methods are considered to reduce the computational time [6].
Whenever possible, it is always recommended to exploit existing symmetries.In particular, the dimension of a problem presenting a helicoidal symmetry, i.e., a combination of translational and rotational symmetries with the same axis, can be reduced from 3D to 2D without loss of accuracy if the calculations are performed in a helicoidal coordinate system.Methods based on this coordinate transformation have first been introduced in optical waveguide simulations [10], [11], [12], and since then applied to electrostatic problems [13], [14], linear magnetodynamic problems [15], [16], [17], and nonlinear magnetodynamic problems with superconducting filaments or tapes [18], [19], [20], [21].
An exact helicoidal symmetry is rarely encountered in practical applications, but different kinds of deformed geometries, curved wires, or conductor organized, e.g., into layers with distinct twist pitch lengths, may exhibit an approximate or partial helicoidal symmetry.Working with helicoidal coordinate systems can still be very useful in such cases, especially in the context of a multi-scale or a sub-problem approach, to compute homogenized parameters that account for the twisting of the filaments (e.g., Rutherford multistrand cables).Furthermore, the 2D helicoidal approach is more accurate than an equivalent 3D approach, as the latter is usually limited in accuracy by non-conformities at element interfaces in unstructured 3D meshes.Extensions and improvements of the helicoidal method are therefore currently being investigated, e.g., in [16], to quantify helicoidal effects in the context of Litz wires.
This paper focuses on the helicoidal transformation method.We start the analysis in Section II by applying the change of coordinates to the h-ϕ-formulation [22], which is an efficient formulation for systems with superconductors [23], and we then state the mathematical conditions for reducing the problem dimension from 3D to 2D.We will refer to the equations resulting from this analysis as the 2D-ξ model, in order to emphasize the fact that it is solved in helicoidal coordinates.As will be shown, a feature of the 2D-ξ model, compared to a conventional 2D model in Cartesian coordinates, is that it solves for fields with three independent components, instead of two.
Depending on the symmetry of the boundary conditions (BC), the study is decomposed in two cases.If the magnetic field excitation is axial and uniform, the BC then also verify the helicoidal symmetry, irrespective of whether there is a transport current or not.The dimension of this problem (geometry plus BC), being helicoidally symmetric, can be reduced from 3D to 2D by simply applying the coordinate transformation method.This approach is not new in the context of superconducting wires [18], [19], [20], but it has not yet been presented with the efficient h-ϕ-formulation.In Section III, the implementation details of this formulation are reviewed with an emphasis on the discretization of a curl-free magnetic field in non-conducting domains.In Section IV, the implementation is verified by comparison with a 3D model in Cartesian coordinates.
If the magnetic field excitation is transverse, then the BC are no longer helicoidally symmetric.A generalization of the method is proposed in Section V for this case.It still results in a 2D model in some situations.To the best of our knowledge, this generalization is a novelty compared to stateof-the-art methods.Attention is again paid to the curl-free property of the magnetic field in non-conducting domains.The generalized model is applied to linear materials in Section VI, and it is shown that it reproduces the predictions of analytical models for the coupling currents [2].We finally provide a brief prospect about the application of the generalized method in the presence of nonlinear materials.
All presented models are implemented in and solved by GetDP [24].Geometry and mesh generation are performed by Gmsh [25].All codes are open-source and available in the Life-HTS toolkit 1 .

II. HELICOIDAL CHANGE OF COORDINATES
Let (x, y, z) be a Cartesian coordinate system.The helicoidal change of coordinates x → ξ and its inverse ξ → x read [11] and respectively, with (ξ 1 , ξ 2 , ξ 3 ) the helicoidal coordinate system.The twisting parameter α ∈ R is the unique parameter of the coordinate transformation, and the pitch length is p = 2π/α.
With this transformation, helices of pitch length p around the z-axis in the Cartesian coordinate system are mapped into straight lines parallel to the ξ 3 -axis in the helicoidal coordinate system.This is illustrated in Fig. 1 with p = 1.A geometry is said to be helicoidally symmetric, or to have a helicoidal symmetry, if there exists a value α for which its description in helicoidal coordinates is ξ 3 -invariant, i.e., independent of ξ 3 .
The Jacobian matrix J of the coordinate transformation Eq. (2) reads with s = sin(αξ 3 ) and c = cos(αξ 3 ).We have det J = 1.The inverse transposed Jacobian matrix J −T , written in terms of the ξ-coordinates, then reads 1 Available: www.life-hts.uliege.be.

A. Helicoidal transformation of fields
The Jacobian matrix describes the mapping of vector components with the transformation.Components of one-forms, such as the magnetic field h, follow the transformation [10], [26] where h x and h ξ denote the components of the field h in the Cartesian and helicoidal coordinate systems, respectively.Components of two-forms, such as the current density j (= curl h), follow the transformation [10], [26] where j x and j ξ denote the components of the field j in the Cartesian and helicoidal coordinate systems, respectively.

B. Problem definition and h-ϕ-formulation
The eddy current problem is governed by the following magnetodynamic (or magneto-quasistatic) equations and constitutive laws [27]: with b, h, j e, µ, and ρ, the magnetic flux density (T), the magnetic field (A/m), the current density (A/m 2 ), the electric field (V/m), the permeability (H/m), and the resistivity (Ωm), respectively.In non-conducting materials, ρ → ∞ and j = 0, and Ampère's law reads In this paper, special attention is paid to satisfy this condition.Type-II irreversible superconductors are characterized by a nonlinear electric response.Assuming isotropy for lowtemperature superconductors, their resistivity is given by the power law [28], where e c = 10 −4 V/m is an electric field threshold defining the critical current density j c (A/m 2 ), and n (-) describes the sharpness of the transition to flux flow.The norm of j is denoted by ∥j∥, with ∥j∥ 2 = j 2 x + j 2 y + j 2 z in the Cartesian coordinates.Finally, all materials are assumed to be nonmagnetic, so that one has µ = µ 0 = 4π × 10 −7 H/m in all domains.
The magnetodynamic problem defined above is solved in a computational domain Ω.Let Ω be a helicoidally symmetric domain.It consists of a conducting domain Ω c made of N connected subdomains, Ω c = ∪ i∈C Ω ci , with C = {1, . . ., N }, surrounded by a non-conducting domain Ω C c .The external boundary of Ω is noted Γ out .Via boundary conditions (BC), the system can be subjected to a given axial magnetic field h axial and/or a given transverse magnetic field h trans .A transport current Īi is imposed to the subdomains Ω ci for i ∈ C I ⊂ C, and a voltage Vi is imposed on the subdomains Ω ci for i ∈ C \ C I = C V .We shall call global conditions (GC) these electric conditions imposed to the conductors of the system.Fig. 2 represents a typical cross section of the problem at hand.We solve the problem defined above with the finite element method.Among the existing finite element formulations, we choose the h-ϕ-formulation [22].It involves the power law written in terms of the resistivity, which has been shown to lead to robust and efficient numerical resolutions for problems involving superconductors characterized by the power law [23].Also, the h-ϕ-formulation strongly verifies the curlfree condition on h in Ω C c , Eq. ( 8), by expressing the magnetic field as the gradient of a scalar potential.This leads to a lower number of degrees of freedom compared to the hformulation [30], that uses instead a spurious non vanishing resistivity to limit the current density in Ω C c .The 3D h-ϕ-formulation reads [23]: from an initial solution at t = 0, find h ∈ H(Ω) such that, for t > 0 and ∀h ′ ∈ H 0 (Ω), we have The integral over Ω of the inner product of f and g is denoted by (f , g) Ω , whereas the operator I i (h) gives the circulation of h around conductor i, which is the net current I i flowing in the conductor.The associated voltage is noted Vi .The function space H(Ω) is the subspace of H(curl; Ω) containing functions that are curl-free in Ω C c and verify the essential BC and the GC [21].The space H 0 (Ω) is the same space as H(Ω) but with homogeneous essential BC and homogeneous GC.For simplicity, we assumed homogeneous natural BC in Eq. (10).

C. The h-ϕ-formulation in helicoidal coordinates
As shown in [31], [32], in order to express the h-ϕformulation Eq. (10) in helicoidal coordinates, it is sufficient to replace the scalar material parameters µ and ρ by the tensors μ and ρ: with the auxiliary tensor T, defined by and its inverse T −1 by This is a consequence of substituting Eqn. ( 5) and (6) into Eq.( 10) and adding a det J factor in the volume integral terms.Beyond these modifications, all calculations can be performed exactly as in Cartesian coordinates [32].
The components of the curl operator in helicoidal coordinates are given by They have the same expression as in Cartesian coordinates, but in terms of the helicoidal coordinates.

D. Conditions for reducing the dimension from 3D to 2D
In the continuous setting, the problems expressed with Cartesian or helicoidal coordinates are equivalent.Indeed, no approximation is introduced and the change of coordinates is regular.For helicoidally symmetric geometries, there are however clear advantages in working with helicoidal coordinates.
First, it involves integrals over domains with ξ 3 -independent sections.
Second, T and T −1 are also ξ 3 -independent, as shown in Eqn.(13) and (14).As a consequence, both the integrand coefficients and the domains of integration in the weak formulation are ξ 3 -independent.
Finally, if the BC on Γ out are also ξ 3 -independent when expressed in helicoidal coordinates, then, the solution h of the h-ϕ-formulation is ξ 3 -independent as well.Hence, the integration along the ξ 3 -direction is trivial and the problem dimension can be reduced from 3D to 2D with a considerable decrease of the computational burden compared to the equivalent 3D problem.
BC are helicoidally symmetric (HS) in the case of a uniform axial magnetic field excitation.In Section III, we describe how the associated 2D problem can be discretized and implemented.We verify the implementation in Section IV.
By contrast, a transverse magnetic field excitation, i.e., a magnetic field in the x-y-plane in the Cartesian coordinate system, does not transform into a ξ 3 -independent field in helicoidal coordinates (see Eq. ( 25)).In this case, the dimension cannot be directly reduced from 3D to 2D.However, simplifications are still possible, eventually also leading to a 2D problem in some situations.We present a novel method for such a situation in Sections V and VI.This method generalizes the case of HS BC, which just becomes a particular case of the general approach.

III. PRACTICAL IMPLEMENTATION OF A FULL
h-ϕ-FORMULATION -HS-BC Starting from Eq. ( 10) with material tensors in Eqn.(11), one could be tempted to implement the h-ϕ-formulation directly as a classical 2D problem with in-plane magnetic field, with the only differences of (i) working in helicoidal coordinates, and (ii) having anisotropic tensors instead of scalar material parameters.But this would not be correct: the fact that the problem is ξ 3 -independent does not imply that the involved magnetic field has only two non-zero (helicoidal or Cartesian) components.
Due to the full anisotropy of tensors μ and ρ, one really has to consider three independent components for the magnetic field h in the h-ϕ-formulation.To emphasize this, we refer to the resulting formulation as a full h-ϕ-formulation in 2D, and we call the associated model the 2D-ξ model.
In this section, we present a practical implementation of this full h-ϕ-formulation.First, we propose a convenient decomposition of the magnetic field, which allows us to reuse the usual function spaces of classical 2D problems.Then, we discuss the discretization of these function spaces.Finally, we explain how to impose the GC and BC.

A. Decomposition of the magnetic field
In the h-ϕ-formulation, the magnetic field h can be decomposed into two parts: an in-plane contribution h ∥ , containing the ξ 1 and ξ 2 -components of h, and an out-of-plane contribution h ⊥ , containing only the ξ 3 -component.We write or, explicitly in terms of their helicoidal components, where h = h(ξ 1 , ξ 2 ) because the solution is ξ 3 -independent.Note that the vectors h ∥ and h ⊥ are not orthogonal.
Because the Jacobian is non-singular, the curl-free condition Eq. ( 8) reads, in the helicoidal coordinate system: from Eq (15) using ∂ ξ3 = 0.With the decomposition defined in Eq. ( 16), the third component of Eq. ( 18) implies that curl h ∥ = 0, which is the same condition as for a classical 2D formulation in which a two-component magnetic field is considered.Then, for the first two components of Eq. ( 18) to be equal to zero, the out-of-plane magnetic field h ⊥ must be uniform in Ω C c .These conditions are introduced in the function space definitions, i.e., they are strongly enforced.They will be made explicit at the space discretization step.
With the explicit decomposition h = h ∥ + h ⊥ , the h-ϕformulation reads as follows.From an initial solution at time where the vectors h ∥ and h ⊥ are coupled by tensors μ and ρ.Note that I i (h ′ ⊥ ) = 0.The function spaces H ∥ (Ω) and H ⊥ (Ω) will be defined in the space discretization step.
For the resistivity in superconducting materials, the power law Eq.( 9) leads to ρ = ρ SC (∥j∥) T. Using Eq. ( 6) and det J = 1, we have, in terms of the components:

B. Space discretization of the magnetic field
Let us consider a finite element mesh for the discretization of the 2D domain Ω, and let us denote by N (Ω i ) and E(Ω i ), the set of nodes and edges, respectively, of the mesh in a given (sub-)domain Ω i , including entities on the boundary of Ω i .
In practice, we can discretize the in-plane magnetic field h ∥ exactly as the two-component magnetic field in a classical 2D h-ϕ-formulation with in-plane magnetic field [23].We use Whitney forms [33]: gradient of node functions w n and cohomology functions c i (cut functions) [34] in Ω C c , and edge functions w e in Ω c \∂Ω c : where coefficients h ∥,e , ϕ n , and I i are the degrees of freedom (DOFs) defining h ∥ in the discrete function space H ∥ (Ω).
We choose to discretize the out-of-plane magnetic field h ⊥ with perpendicular edge functions w n = w n êξ3 , associated with nodes.To account for the fact that h ⊥ must be uniform in each region of Ω C c , we introduce global functions in Ω C c .Let K be the number of connected regions in Ω C c .We describe the out-of-plane magnetic field with the expansion with w n the perpendicular edge function associated with node n in N (Ω c \ ∂Ω c ), and p i a global shape function defined as the sum of all perpendicular edge functions associated with nodes in the i th connected region of Ω C c , including those on its boundary, for i ∈ {1, . . ., K}.The support of the shape function p i is therefore not restricted to Ω C c : it is non-zero on a layer of one element adjacent to ∂Ω c in Ω c .This defines the discrete function space H ⊥ (Ω), with DOFs h ⊥,n and D i .Both h ∥ and h ⊥ are described by discrete 1-forms, and so is their sum, h.
For simplicity, in the following, we assume that there is only one connected non-conducting region Ω C c , the exterior of the wire, such that K = 1, and we rename D 1 = D.

C. Global conditions and boundary conditions
For the GC, a current Īi , for i ∈ C I , can be imposed exactly as in a classical 2D h-ϕ-formulation with in-plane magnetic field, i.e., strongly via the degree of freedom I i associated with the cut function c i for the corresponding conducting domain Ω ci .Alternatively, an applied voltage Vi , for i ∈ C V , can be imposed weakly in the global term of the formulation Eq. (19).
For the BC, we consider a circular external boundary Γ out , placed in Ω C c sufficiently far from the conductors such that we can assume ∂ t b • n| Γout = 0, with n the outer normal vector.This condition is implicitly imposed for h ∥ in Eq. ( 19) with homogeneous natural BC on Γ out .This lets the z-component of the magnetic field, h z , undetermined on Γ out .It corresponds to the axial magnetic field, that we can freely impose.We derive below how to translate this into a BC on h ⊥ | Γout in helicoidal coordinates.
Let us first consider the situation with a zero axial magnetic field.At a sufficiently large distance R out from the center of conductors carrying a total net current intensity I, the magnetic field tends to be purely azimuthal and axisymmetric.We have h x = I 2πRout (− sin θ cos θ 0) T , with θ = atan2(y, x).In terms of the helicoidal coordinates, on the plane ξ 3 = 0, it reads using ξ 2 = R out sin θ and ξ 1 = R out cos θ for ξ 3 = 0. Consequently, to satisfy h z | Γout = 0, one has to impose that h ξ3 | Γout = Iα/2π.This can be done by fixing the degree of freedom D associated with the basis function p in Ω C c in Eq. ( 22) to the value D = Iα/2π.Note that this value does not depend on R out .
By superposition, if one wants to impose a non-zero axial magnetic field h axial on the external boundary Γ out in addition to a net current intensity I, we can impose the following condition: because the axial magnetic field Cartesian components h x = (0 0 h axial ) T transform into h ξ = J T h x = (0 0 h axial ) T in helicoidal coordinates.

IV. VERIFICATION AND APPLICATION -HS-BC
In this section, we first compare the solution of the 2D-ξ model in helicoidal coordinates to the solution of a classical 3D h-ϕ-formulation on a simple problem in order to verify the implementation.We also quantify the computational gain offered by reducing the dimension from 3D to 2D.Then, we apply the 2D-ξ model on a more involved geometry to illustrate the capabilities of the approach.

A. Verification problem
We consider a wire made of six identical Nb-Ti superconducting filaments, twisted and embedded in a copper (Cu) matrix, as illustrated in Fig. 3.In order to simplify the geometry, the cross sections of the filaments are assumed to be disks, and the 3D geometry is generated by a helicoidal extrusion of them.This is of course an approximation of a realistic geometry.If needed, cross sections of round twisted filaments can be computed accurately using envelope theory as in [15] or CAD tools [25] as in [17].We assume that Nb-Ti resistivity is characterized by Eq. ( 9) with constant and uniform j c = 7×10 9 A/m 2 and n = 50, and that the copper resistivity is ρ Cu = 1.81 × 10 −10 Ωm.There is no insulation between the filaments and the matrix, so that the wire behaves as a single conducting cylinder.A net transport current I(t) = 0.5 I c sin(2πt/T ) is imposed in the wire, with T = 0.1 s and I c = 162 A, and we impose h axial = 0 A/m.

B. Implementation of the 3D model
We consider the 3D geometry represented in Fig. 3(b).It represents a periodic cell of one-sixth of a whole pitch length p.Note that building and meshing the 3D model represented in Fig. 3(b) is not a trivial task.To account for the periodicity of the problem, the mesh must be identical on the top and bottom boundaries of the domain; as an h-ϕ-formulation is used, cohomology basis functions must also be periodic.The quality of the mesh inside the filaments plays an important role for the accuracy of the resulting numerical solution.We observed that better results are obtained with a structured mesh inside the filaments.Generating the mesh with such constraints is possible with Gmsh [25].The periodic support for the cohomology basis function is generated as described in [34], [35] and illustrated in Fig. 4. We set a homogeneous natural BC on the external boundary Γ out , so that ∂ t b • n| Γout = 0 is weakly enforced.For the top and bottom boundaries Γ up and Γ down , which are topologically identical, the periodic condition h × n| Γup = −h × n| Γdown is imposed.On conducting boundaries ∂Ω c ∩(Γ up ∪Γ down ), this is done by forcing the equality of the degrees of freedom associated with topologically identical edges of these boundaries.On non-conducting boundaries ∂Ω C c ∩ (Γ up ∪ Γ down ), the periodic constraint is enforced via the magnetic scalar potential.We impose ϕ| Γup = ϕ| Γdown + h axial p/6.The total current intensity flowing in the conducting domain made up of the filaments and matrix is imposed via the (periodic) cohomology basis function whose generating edges are highlighted in Fig. 4.
Note that in the present case of HS-BC (transport current or axial field), the 3D reference model could be defined on a length shorter than p/6 along z, if one adapts the periodic mesh and the periodic cut accordingly.We chose a length of p/6 so that the reference model will also be valid in the transverse field case.
Before comparing the results, we first verify that the 3D model indeed produces a helicoidally symmetric solution.For illustration, from the 3D numerical solution, we extract the magnetic field h and the current density j along the helicoidal fiber of pitch length p passing at point x = a, b, 0 , with a = 180 mm and b = 11 mm, from z = 0 to z = p (see Fig. 4).We exploit the periodicity of the problem to obtain values for z > p/6.The Cartesian and helicoidal components of vectors h and j are represented in Fig. 5 for a relatively fine tetrahedral mesh (144 870 DOFs), at time t = T /4.Helicoidal components are obtained using the one and two-forms transformation relations, Eqn. ( 5) and ( 6).The oscillations and spikes along the fiber represent interelement non-conformities, which are expected with lowestorder tetrahedral Whitney shape functions.These oscillations decrease in amplitude with mesh refinement.Up to these inter-element variations, the 3D solution correctly presents a helicoidal symmetry.It is also interesting to notice that the current density has non-zero ξ 1 and ξ 2 -components, and that the ξ 3 -component of the magnetic field is not equal to zero.This illustrates the need for a three-component magnetic field in the 2D helicoidal model.

C. Comparison of the results from the 3D and 2D-ξ models
We now compare the results of the 2D-ξ problem in helicoidal coordinates with the reference 3D problem described above.Note that for the 2D-ξ model, in this particular case, we could further exploit the symmetry and model only one-sixth of the circular region depicted in Fig. 3(a) using periodic BC on the symmetry boundaries as well as an adapted cohomology function in Ω C c , hence reducing the computational cost even more.We however choose to model the full 2D cross section.
The solution of the 2D-ξ model on a medium mesh resolution (4 700 DOFs) is represented in Fig. 6.The current mostly flows in the superconducting filaments, as shown by the different scales for the middle and right subfigures.On the left subfigure, one can see that the current flow in the twisted filaments induces a non-zero z-component h z of the magnetic field at the center of the wire.A comparison of the local magnetic field of the 2D-ξ model with that of the reference 3D model is given in Fig. 7, along the dashed red line highlighted in Fig. 6, for two mesh resolutions.The solution of the 3D model is taken on the plane z = ξ 3 = 0, but this choice is arbitrary: as was shown in Fig. 5, up to the inter-element variation, the solution of the 3D model is also ξ 3 -independent.Solutions of the 2D-ξ and 3D models match locally.We verified and this is also the case for the current density (not represented in the figures).A comparison of the AC loss is given in Fig. 8.The AC loss per unit length along êz in both the superconducting filaments and the conducting matrix are compared for the two models, and for two mesh resolutions.For the 2D-ξ model in helicoidal coordinates, the AC loss is computed as ρ j ξ , j ξ Ωc , where Ω c is either restricted to the filaments, or to the matrix.For the 3D model, the integral (ρ j x , j x ) Ωc is computed over the 3D domain with Cartesian coordinate system, and the result is divided by p/6, to obtain the AC loss per unit length as well.Note that both models include all loss contributions by construction: hysteresis losses in the filaments, as well as coupling and eddy current losses in the matrix [36].Meshes for the coarse resolution in the z = ξ 3 = 0 plane are similar for the 2D-ξ and 3D models, as well as meshes for the fine resolution.However, we observed that the 2D-ξ model solution is less sensitive to the mesh resolution.This is due to the inter-element non-conformities in a tetrahedral 3D mesh, which should be made significantly lower (by mesh refinement) to ensure an accurate evaluation of the quadratic quantity representing the AC losses.Meshes with prisms, i.e., extruded triangles, in the filaments were also tested.They give slightly better results, but also increase the complexity of the meshing step, as pyramids must be used as transition elements between prisms in the filaments and tetrahedra outside of them.
The local and global quantity agreement shows the validity of the 2D-ξ model in helicoidal coordinates.The dimension reduction allows for a very large reduction of the computational cost.This is demonstrated in Table I, which compares the performance of the 2D-ξ and 3D models on meshes with similar characteristic length for the finite elements (triangles in 2D-ξ and tetrahedra in 3D).The fine 2D-ξ model is more than two orders of magnitude faster to solve than the fine 3D model.

D. Application to a 54-filament wire
As a more realistic geometry, we consider a wire with 54 filaments arranged in a hexagonal lattice with filament center spacing of d = 110 µm, as represented in Fig. 9    Note that during the first transport current increase (for t < T /4), no filament carries a negative current.This is in contradiction to what was obtained in [37] with an alternative method on a similar problem, where a critical state model is considered and the current density is assumed localized along a line within every filament, which therefore leads to neglecting the spatial extension of the filaments.

V. EXTENSION TO NON-HELICOIDALLY-INVARIANT BOUNDARY CONDITIONS -GENERAL BC
When BC are not HS, the dimension of the problem cannot be directly reduced from 3D to 2D on basis of the geometrical symmetry only.This is the case when a wire is subjected to a uniform transverse magnetic field.For an applied magnetic field h x = (0 1 0) T , we have which is not ξ 3 -independent, see Fig. 12.As a consequence, the solution of the magnetodynamic problem will not be ξ 3independent either.The periodic structure of the problem can however be exploited by expressing the solution as a series of periodic functions with respect to ξ 3 .
In this section, we present this approach and show that it generalizes the method described in Section III.In particular, we show that it also leads to a 2D model in helicoidal coordinates, which has the potential of considerably reducing the computational cost compared to a 3D model.

A. Fourier decomposition of the magnetic field
Given the p-periodicity with respect to ξ 3 , by separation of variables, we can expand the magnetic field h = h(ξ 1 , ξ 2 , ξ 3 ) in the following series: with the modes f k = f k (ξ 3 ) that are functions of ξ 3 only and that are defined as and with the spatial Fourier coefficients h k = h k (ξ 1 , ξ 2 ) that are three-component vector functions of ξ 1 and ξ 2 .
The modes f k are mutually orthogonal and have a unit norm, denoted as ∥f k ∥ = 1, in the sense of the following inner product: They also satisfy the following property: Introducing a decomposition of the magnetic field into its in-plane and out-of-plane components as in the HS case of Section III for each Fourier coefficient h k , we can rewrite Eq. ( 26) as with the h ∥,k containing the ξ 1 and ξ 2 -components of h k , and the h ⊥,k containing its ξ 3 -component.Equation ( 30) actually generalizes the decomposition in Eq. ( 16) for the case of HS BC.Indeed, in the case of HS BC, the only mode that is involved is f 0 (ξ 3 ) = 1, with coefficients h ∥,0 = h ∥ and h ⊥,0 = h ⊥ , and the coefficients of the other modes, h ∥,k and h ⊥,k , ∀k ∈ Z 0 , are all equal to zero.

B. Space discretization with curl-free functions in Ω C c
The curl of decomposition (30) reads where êξ3 is the unit vector in the ξ 3 -direction.The only term in Eq. ( 31) contributing to the ξ 3 -component of the curl involves the curl of h ∥,k .We can therefore keep the same discrete function space for the h ∥,k as for h ∥ in Section III, however without the i I i c i term of Eq. ( 21) for k ̸ = 0, as transport currents only contribute to the fundamental mode with f 0 (ξ 3 ) = 1.
As in the ξ 3 -independent case, we express the out-of-plane magnetic field h ⊥,k as a sum of perpendicular edge functions.
But now, the h ∥,k functions also contribute to the ξ 1 and ξ 2components of the curl of h for k ̸ = 0 in Eq. ( 31) via the cross product term.Therefore, the curl-free condition in Ω C c is no longer met with a uniform out-of-plane magnetic field in Ω C c , for k ̸ = 0. Instead, as is shown below, the curl-free condition induces a coupling between the in-plane and outof-plane magnetic field contributions in Ω C c .For simplicity, as was done before, we assume that there is only one connected non-conducting region Ω C c .Using curl-free in-plane functions h ∥,k in Ω C c and Eq. ( 31), the curl-free condition on Using the mode property Eq. ( 29), this yields which results in the following condition, ∀k ∈ Z: For k = 0, we retrieve the same condition as in the helicoidally symmetric problem, that is, h ⊥,0 must be uniform in Ω C c , with a value given by Eq. ( 24).For k ̸ = 0, the condition can be enforced via the independent degrees of freedom of the inplane and out-of-plane magnetic field contributions.Indeed, in Ω C c , we have the expansions where w n êξ3 = w n is the perpendicular edge function of node n, with w n the usual node function.In terms of the individual degrees of freedom, Eq. ( 34) reads This equation is valid over the whole domain Ω C c if and only if the first parenthesis is constant.This is the case if, for k ̸ = 0, That is, to ensure a curl-free magnetic field in Ω C c , the degrees of freedom of the mode h ⊥,k must be linked directly to those of the mode h ∥,−k in Ω C c (or vice-versa).This link between the degrees of freedom strongly ensures that curl h = 0 in Ω C c and allows for a significant reduction of the number of unknowns, hence a reduction of the computational cost of the resolutions.

C. Boundary conditions for a transverse magnetic field
The transverse magnetic field defined in Eq. ( 25) applied as a BC on Γ out only involves the modes f −1 (ξ 3 ) and f 1 (ξ 3 ).We have, in helicoidal components: We can verify that they satisfy Eq. ( 34).

D. Derivation of the h-ϕ-formulation with linear materials
With linear materials, orthogonality allows solving modes with different values of |k| (i.e., including −k and k) independently.For each value of |k|, the integration along ξ 3 gives an independent set of equations, written in terms of the unknown Fourier coefficients h ∥,−k , h ⊥,−k , h ∥,k , and h ⊥,k .These coefficients are functions of ξ 1 and ξ 2 only, and hence the problem is 2D.The formulation is derived in the Appendix.
For nonlinear materials, the modes are no longer decoupled.We provide observations and comments on how to handle this situation in Section VI-C.

VI. VERIFICATION AND APPLICATION -GENERAL BC
In this section, we first verify the implementation of the generalized 2D-ξ method by comparing its results with those of a 3D reference model, for linear materials.We then apply the method on a 54-filament wire and discuss the different contributions to the total AC loss, still with linear materials.Finally, we comment on the application of the method in the case of nonlinear materials.

A. Verification with linear materials
The validity of the approach with linear materials is verified by comparing the results of the 2D-ξ model with those obtained with a classical 3D model.We consider the same geometry as in Section IV, but with a constant resistivity in the filaments, and with a uniform transverse magnetic field instead of an imposed transport current.
The filaments have a constant resistivity ρ SC = 3.3 × 10 −14 Ωm (dummy value chosen for verification), and the matrix has a constant resistivity ρ Cu = 1.81 × 10 −10 Ωm.The system is subjected to a transverse magnetic field along y, increasing from 0 T to 0.1 T with a constant ramp-up rate of 18 T/s.
Comparisons with the solution of the 3D problem are given in Figs. 14 and 15, along a characteristic line in the z = ξ 3 = 0 plane and along a helicoidal fiber of pitch length p, passing at point x = r, 0, 0 , with r = R ℓ + 0.8R f , from z = 0 to z = p.Both models agree with each other.
Fig. 14: Magnetic field components along the dashed line represented in Fig. 13, at z = 0, for the 3D and 2D-ξ models with a fine mesh resolution for µ0hy = 0.1 T.
As in the HS-BC case, exploiting the geometrical symmetry allows for a strong reduction of the computational work.It should however be mentioned that the 2D-ξ model with transverse field BC involves double number of degrees of freedom compared to the same model with HS-BC, as two modes are needed to represent the transverse field (−k and k, compared to k = 0 only).
The 2D-ξ model still leads to a considerable reduction of DOFs compared to the 3D model.Indeed, taking values of the fine mesh resolution from Table .I, the 3D model involves 145k DOFs whereas the 2D-ξ model with two modes only involves 12k DOFs.0 0.2 0.4 0.6 0.8 1 z (mm) Fig. 15: Magnetic field along the helicoidal fiber of pitch length p, passing at point x = r, 0, 0 , with r = R ℓ + 0.8Rf (represented by the red dot in Fig. 13), from z = 0 to z = p, for the 3D and 2D-ξ models for µ0hy = 0.1 T. (Left) Cartesian components of the vectors.(Right) Helicoidal components of the vectors.

B. Application on a 54-filament wire with linear materials
We now consider the 54-filament geometry defined in Fig. 9 but with a linear material in the filaments.We fix the resistivity in the filaments to ρ SC = 1.81 × 10 −15 Ωm and in the matrix to ρ Cu = 1.81 × 10 −10 Ωm.
Choosing such a low resistivity in filaments leads to an approximated model for superconducting wires at low magnetic fields (below filament saturation).We will see that this linear model reproduces the coupling current dynamics observed in the copper matrix of superconducting wires.The validity of the linear model is however limited to this.It does not describe superconducting hysteresis effects in filaments, and hence does not allow for superconducting loss calculation.Instead, in the following, the computed losses in the filament region will be those of a normal resistive material.
We compute the AC loss induced by a time-varying transverse magnetic field µ 0 h = b max sin(ωt)ê y , with b max = 0.1 T, as a function of the frequency f = ω/2π.As in the previous section, BC are such that only the modes f −1 (ξ 3 ) and f +1 (ξ 3 ) are excited (see Eqn. (39) to (42)).
Moreover, because the materials are linear and the excitation is harmonic, the problem can be solved in the frequency domain.To this end, we write the problem in terms of the auxiliary complex quantity ĥ(ξ), the phasor of the magnetic field.The phasor is related to the physical magnetic field by h(ξ, t) = ℜ ĥ(ξ)e iωt , with i = √ −1, and we replace all time derivatives in the formulation by a multiplication by iω.
The time-average instantaneous loss density, in W/m 3 , reads, in terms of Cartesian and helicoidal components of the phasor ĵ for the current density: where ĵ⋆ denotes the tranposed complex conjugate of ĵ.Note that both sides of Eq. ( 44) are real since ρ is a scalar and ρ is a hermitian tensor.The total loss per unit length is obtained by integrating Eq. ( 44) over the whole wire cross section.
We decompose this total loss into separate contributions, which allows for an easier interpretation of the results.The filament loss is the integral of Eq. ( 44) on the filament region only.The coupling loss is the integral of Eq. (44) on the matrix region only, taking only the in-plane components of ĵx into account (the x and y Cartesian components).Finally, the eddy current loss is the same but with only the out-of-plane component ĵz .
We present the results for frequencies ranging from 10 −2 Hz to 10 5 Hz in Fig. 16 for a pitch length p = 10 mm.Values from a 3D reference model are also given for comparison, the agreement with the 2D-ξ model is very good.The current distribution at two distinct frequencies is shown in Fig. 17     Figure 16 shows that the dominant loss contribution depends on the frequency.This can be interpreted as follows.
Neglecting the effect of the twist, we expect the peak of the filament loss to arise when the diffusion skin depth δ SC = 2ρ SC /ωµ is comparable with the radius of the wire R f .We have δ SC /R f = 1 for the frequency f = 0.22 Hz, which is not too far from the peak in Fig. 16(b).
A change of regime for the eddy current loss per cycle should arise when the skin effect starts to play a role in the matrix.This is expected to happen when the diffusion skin depth δ Cu = 2ρ Cu /ωµ is comparable with the thickness of the outer sheath of the matrix D os ≈ 80 µm.Here, we have δ Cu /D os = 1 for the frequency f = 7.2 kHz, which coincides with the peak in Fig. 16(b).Below this frequency, in the 0.1 kHz to 5 kHz range, most of the field is shielded by currents in the filaments, that are coupled via coupling currents, as discussed below.
Coupling losses are due to currents flowing between the filaments, known as the coupling currents [38].They are represented by the arrows in Fig. 17.It is worth mentioning that they are on average flowing anti-parallel to the applied magnetic field, as predicted by analytical models [38], [36].Their dynamics is that of an RL-circuit governed by a time constant τ c and they contribute to a loss per cycle and per unit length q cycle (J/m).Simplified models propose [36]: with ρ eff the effective resistivity of the matrix, accounting for the presence of the filaments [2].In the present case in which we assume no insulation between the filaments and the matrix, assuming the filaments have negligible resistivity, we can estimate ρ eff as follows [39]: with λ the filling factor of the filaments in the wire.Here, λ = 0.44 so that τ c = 23 ms.The associated frequency is f c = (2πτ c ) −1 = 7 Hz, which roughly corresponds to the position of the peak value of the coupling loss per cycle in Fig. 16(b).Below the peak frequency, the filaments are mostly decoupled, as the magnetic field does not change fast enough for large coupling currents to appear.Above the peak frequency, they get more and more coupled, as illustrated in Fig. 18.
As an illustration of the effect of p on the coupling losses, we give in Fig. 19 the total and coupling losses for different values of the pitch length.We can verify the agreement with the analytical prediction Eq. ( 45): the peak position of the coupling losses scales quadratically with p, affecting the total loss significantly.The curve for pure copper cylindrical conductor of the same radius R w , with ρ SC = ρ Cu = 1.81 × 10 −10 Ωm, is given for comparison.
Note that for the pure copper case, we have δ Cu /R w = 1 for f = 183 Hz, which roughly corresponds to the position of the peak of AC loss per cycle.Frequency (Hz) AC loss per cycle (J/m) Fig. 19: AC loss as a function of the frequency of an applied transverse magnetic field for linear materials and different pitch lengths.Solid curves represent the total loss.Dashed curves represent the coupling loss only.The main effect of the twist is to shift the coupling loss curves to higher frequencies for decreasing values of p.The total loss for a wire made of copper only (2D model) is given for comparison (dotted curve).The legend is valid for both subfigures.
Figure 19 clearly shows the beneficial effect of twisting the filaments at low frequencies.Wires with smaller pitch length indeed have shorter time constants and are less subject to coupling current losses at low frequencies.It must be mentionned that the twist however does not reduce loss for all frequencies, which is in agreement with experimental measurements, e.g., in [40].
The simple linear model discussed here allows for a qualitative description of coupling current losses that are representative of real superconducting wires for low applied magnetic field only.With superconducting filaments, saturation effects change the coupling current dynamics for higher field amplitudes [2].Such effects cannot be reproduced with a linear model.
As already said, the analysis of this linear model must therefore be carried out with caution.Total loss evaluations for superconducting wires cannot be extracted from this model as hysteresis losses of superconducting filaments are not part of the linear model.The inclusion of nonlinear material properties is necessary for such analysis.In the next section, we present the challenges that such an inclusion brings to the helicoidal transformation approach with general BC.

C. Comments for nonlinear materials
In the presence of nonlinear materials, such as superconducting filaments with a power law resistivity described by Eq. ( 9), mode decoupling is no longer possible general BC.As derived in the Appendix, the eddy current term of the formulation expands as a double sum on k, k ′ ∈ Z of the terms given by Eq. (51).And each term in Eq. (51) involves the tensor ρ.For a superconducting filament, this tensor depends on the full local current density, which couples the modes with different values of |k|.A large number of modes in Eq. ( 30) is therefore likely to be excited by a transverse magnetic field.
As the integral along the ξ 3 -direction can no longer be computed a priori, the resulting problem is no longer twodimensional, which makes it qualitatively different from the 2D-ξ model with linear materials.
To assess the importance of this mode coupling, we can use the 3D model.We show in Fig. 20 the evolution of the magnetic field and the current density along one helicoidal fiber, obtained with the 3D model with the same material parameters as in Section IV, but subject to a transverse magnetic field.As can be seen on the bottom-right plot, the magnetic field in helicoidal coordinates cannot be described only with the two modes f −1 (ξ 3 ) and f +1 (ξ 3 ) as in the linear case.Higher modes are excited.
The amplitude of the different modes can be quantified by a discrete Fourier transform of the magnetic field evolution along this helicoidal fiber.This is illustrated in Fig. 21.Small, but non-negligible contributions are brought by modes |k| > 1.
Whether a description with a limited number of modes would lead to satisfying evaluations of losses or not is not an obvious question; knowing a priori how many modes should be considered on a new geometry remains an open question.Further investigations in that direction are necessary.

VII. CONCLUSION
In this work, we applied a change of coordinates on the h-ϕ-formulation for modelling multifilamentary wires presenting a helicoidal symmetry.This led to a reduction of the geometrical dimension from 3D to 2D, hence allowing for a substantial gain in terms of computational effort.We separated the study in two steps, depending on the helicoidal symmetry of the boundary conditions (BC).In both cases, we described in details the spatial discretization of finite element fields in helicoidal coordinates.In particular, we emphasized the necessity of using three independent components for the unknown fields.We then successfully verified our implementation against standard 3D models.
In the case with no external field (e.g.transport current situation only) or with an axial magnetic field, BC are helicoidally symmetric and the method can be directly applied to nonlinear materials.The approach is exact in the sense that no approximation is introduced in the continuous setting.The proposed method can be directly applied on single-layer CORC ® cables [41] or twisted stacked-tape conductors [42].
In the case of a transverse magnetic field, BC are no longer helicoidally symmetric, but the approach was generalized and applied to linear materials.We presented a study of coupling current induced losses in the harmonic regime, and we finally commented on a possible extension to nonlinear materials.For nonlinear materials with general BC, further investigations are necessary.
In the linear case in which ρ is not a function of the fields, we can integrate each term along the geometry invariant ξ 3direction over one pitch length p, divide by p, use the mode property Eq. ( 29), and exploit the mode orthonormality.For k = 0, because d ξ3 f 0 = 0, only terms for k ′ = 0 survive, and they are decoupled from all other terms (k ̸ = 0).These terms are the same as the ones implemented in the case of helicoidally symmetric boundary conditions (HS-BC): ρ curl h ∥,0 , curl h ′ ∥,0 Ωc + ρ curl h ⊥,0 , curl h ′ ∥,0 Ωc + ρ curl h ∥,0 , curl h ′ ⊥,0 Ωc + ρ curl h ⊥,0 , curl h ′ ⊥,0 Ωc .
For k ̸ = 0, only one term of the sum on k ′ survives for each term, either k ′ = k, or k ′ = −k.Indeed, Eq. ( 29) induces the coupling of the modes k and −k.For a given value of k ̸ = 0, in Eq. (51), the only terms that remain are To these terms, another set needs to be added, with the opposite value of k, k ⋆ = −k.In total, this gives eighteen individual terms for the eddy current contribution, for each value of |k| ̸ = 0.

Fig. 3 :
Fig. 3: Wire geometry for the verification of the helicoidal transformation consisting of six twisted Nb-Ti filaments embedded in a copper matrix.(a) Geometry in a ξ1-ξ2 plane (or in the x-y plane for z = 0).(b) One-sixth of a pitch length.The copper matrix in (b) is not represented, for clarity.The filaments have a radius of R f = 35 µm and their centers are at a distance R ℓ = 98 µm from the center of the wire.The wire has a radius of R w = 155 µm and a pitch length of p = 1 mm.The air is modelled outside of the wire up to a distance R out = 500 µm.We assume that Nb-Ti resistivity is characterized by Eq. (9) with constant and uniform j c = 7×10 9 A/m 2 and n = 50, and that the copper resistivity is ρ Cu = 1.81 × 10 −10 Ωm.There is no insulation between the filaments and the matrix, so that the wire behaves as a single conducting cylinder.A net transport current I(t) = 0.5 I c sin(2πt/T ) is imposed in the wire, with T = 0.1 s and I c = 162 A, and we impose h axial = 0 A/m.

Fig. 4 :
Fig. 4: Periodic support for the cohomology basis function to impose a transport current I(t) in the 3D verification model.The red curve in the filaments is a portion of the helicoidal fiber along which the solution is represented in Fig. 5.

Fig. 6 :
Fig. 6: Magnetic field (left), current density in the filaments (middle), and current density in the matrix (right) at time t = T /4 for the 2Dξ problem solved with the helicoidal coordinate system.The arrows represent the in-plane x and y-components of h and j, whereas the triangular elements are colored as a function of the out-of-plane zcomponent of h and j.The dashed red line in the left figure is the cut along which the magnetic field is represented in Fig. 7.

Fig. 7 :
Fig.7: Magnetic field along the dashed red line represented in Fig.6, for the 3D and 2D-ξ models, at time t = T /4, with coarse (up) and fine (down) mesh resolutions.

Fig. 8 :
Fig.8: AC losses in the superconducting filaments (up) and in the conducting matrix (down) for a transport current I(t), as a function of time, for two mesh resolutions, with the 2D-ξ model in helicoidal coordinates and the 3D verification model.

Fig. 10
Fig.10shows the time evolution of the current in the filaments depending on their position.As expected, the current density progressively penetrates into inner layers of the wire.Due to the twist, the current flowing in the outer filaments generates a non-zero h z component inside the wire.Circulating in-plane currents therefore appear in the inner layers to shield this axial magnetic field, as illustrated in Fig.11.

Fig. 10 :
Fig. 10: Distribution of the current among the filaments, as a function of their distance to the center O, for the 54-filament geometry.

Fig. 11 :
Fig. 11: Magnetic field (left) and current density (right) for the 54filament wire at time t = 0.1 s.Only one quarter of the geometry is shown.Arrows represent the in-plane components and the elements are colored as a function of the value of the z-component of the vectors.Note that two different scales are used for the magnetic field for clarity.

Fig. 12 :
Fig.12: Uniform magnetic field along êy for the six-filament geometry represented in Fig.3.(Left) In the physical space on the plane x = 0. (Right) In the helicoidal coordinate system represented as an orthogonal system on the plane ξ1 = 0.

Fig. 13 :
Fig. 13: Solution of the 2D-ξ model with linear materials, on the z = 0 plane, for a transverse field of 0.1 T. The arrows represent µ0h, and the triangular elements are colored as a function of the value of µ0hz, using the color map on the top.The dashed red line is where the field is taken for Fig. 14, and the red dot along that line represents the intersection with the plane z = 0 of the helicoidal fiber along which the field is taken for Fig. 15.(Left) Mode with f−1(ξ3) (see Eq. (43)).(Right) Mode with f+1(ξ3) (see Eq. (43)). .

Fig. 16 :
Fig. 16: AC loss as a function of the frequency of an external transverse magnetic field for linear materials and p = 10 mm.The total loss as well as separate contributions (filament, coupling, eddy) are shown.The legend is valid for both subfigures.The markers denote the total loss obtained by a 3D model, for verification.

7 Fig. 17 :
Fig. 17: Real part of the current density distribution in the matrix in transverse magnetic field in harmonic regime.Arrows represent the in-plane x-y-components, and elements are colored as a function of the out-of-plane z-component.

Fig. 18 :
Fig. 18: Real part of the current distribution in the filaments in transverse magnetic in harmonic regime.Arrows represent the inplane x-y-components, and elements are colored as a function of the out-of-plane z-component.

µ0h ξ 2 µ0h ξ 3 Fig. 20 :hξ 1 µ0 hξ 2 µ0 hξ 3 Fig. 21 :
Fig.20: Current density (up) and magnetic field (down) along the helicoidal fiber of pitch length p, passing through point x = r, 0, 0 , with r = R ℓ + 0.8Rf from z = 0 to p, for a transverse applied magnetic field along êy and Nb-Ti filaments.(Left) Three components of the vectors in the x-space.(Right) Three components of the vectors in the ξ-space.Solution of the 3D model on a fine mesh with prismatic element in the filaments.

TABLE I :
. Filament radius is R f = 45 µm, wire radius is R w = 500 µm, and Performance comparison for the 3D and 2D-ξ models with imposed current and no axial magnetic field, computed with 150 time steps from t = 0 to t = 5 T /4 on a single Intel Core i7 2.2 GHz CPU.DOF: degree of freedom.It.: iteration (Newton-Raphson).