2D Fourier finite element formulation for magnetostatics in curvilinear coordinates with a symmetry direction
22D Fourier finite element formulation for magnetostaticsin curvilinear coordinates with a symmetry direction
Christopher G. Albert , Oszkár Bíró , Patrick Lainer Max-Planck-Institut für Plasmaphysik, Boltzmannstraße 2, 85748 Garching, GermanyGraz University of Technology – Graz Center of Computational Engineering Institute of Theoretical and Computational Physics, Petersgasse 16, 8010 Graz, Austria Institute of Fundamentals and Theory in Electrical Engineering, Inffeldgasse 18, 8010 Graz, Austria
Abstract
We present a numerical method for the solution of linear magnetostatic problems in domains with a symmetrydirection, including axial and translational symmetry. The approach uses a Fourier series decomposition of thevector potential formulation along the symmetry direction and covers both, zeroth (non-oscillatory) and non-zero(oscillatory) harmonics. For the latter it is possible to eliminate one component of the vector potential resultingin a fully transverse vector potential orthogonal to the transverse magnetic field. In addition to the Poisson-likeequation for the longitudinal component of the non-oscillatory problem, a general curl-curl Helmholtz equationresults for the transverse problem covering both, non-oscillatory and oscillatory case. The derivation is performedin the covariant formalism for curvilinear coordinates with a tensorial permeability and symmetry restrictions onmetric and permeability tensor. The resulting variational forms are treated by the usual nodal finite element methodfor the longitudinal problem and by a two-dimensional edge element method for the transverse problem. Thenumerical solution can be computed independently for each harmonic which is favourable with regard to memoryusage and parallelisation.
I. I
NTRODUCTION
Two-dimensional formulations of electrodynamic problems in translationally and axisymmetric systemsare commonly used for numerical treatment by the finite element method [1, 2]. This includes simplifiedmodels of three-dimensional systems, as the number of degrees of freedom is substantially reducedcompared to the full model. In the simplest case all quantities are assumed to be independent of thesymmetry variable, thereby reducing the problem dimension by one. In the more general Fourier finiteelement method also oscillatory field components are treated by combination of a Fourier expansion inthe symmetry coordinate with a numerical solution for each individual harmonic (see, e.g., [3, 4, 5] forthe treatment of Maxwell’s equations). For linear problems the Fourier finite element method has theadvantage of trivial parallelisability over individual harmonics in comparison to a 3D computation.Here we treat the magnetostatic problem using the vector potential curl - curl formulation in domainswith a symmetry direction, including in particular the translationally and axisymmetric case. Due to thelinearity of the curl operator, it is possible to choose a different gauge for each Fourier harmonic andobtain a valid vector potential as the total sum. A similar approach has been taken for the divergence-freeinterpolation of a given magnetic field in [6]. As a result, the oscillatory part of the curl - curl equationcan be treated as a a fully two-dimensional problem using a transverse vector potential. The resultingequations can be viewed as a generalisation of the transverse part of the axisymmetric problem and are in Preprint version: September 1, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] A ug analogy to the transverse formulation of waveguide models [7]. The approach is implemented using themetric-free covariant formulation of Maxwell’s equations in curvilinear coordinates (see [8, 9] for detailsor [10] for a concise introduction) that is briefly restated in the notation used for further derivations. Herelinear constitutive relations contain the permeability tensor represented by its covariant density componentsimplicitly including the influence from curvilinear coordinates via their respective metric tensor. Whilebeing mathematically equivalent to the formulation via differential forms the covariant notation remainscloser to the traditional formulation via differential operators in flat space. For convenience, Table II liststhe relevant designations for both approaches.Explicit expressions with a scalar permeability are given for the Cartesian case with π -periodicity in z -direction and axisymmetric systems in cylindrical coordinates. Application to other axisymmetric systemssuch as spherical, toroidal or symmetry flux coordinates [11] is straightforward by inserting the componentsof their respective metric tensor. As long as the latter fulfils the necessary requirements, even moreexotic coordinate systems (e.g. helical [12]) should be realisable. The resulting variational problems arediscretised in coordinate space using a two-dimensional (Nédélec) edge finite element formulation [13, 14]of lowest order for the transverse equations and first order Lagrange elements for the longitudinal non-oscillatory Poisson-like equation. A convergence study for benchmarking problems with existing analyticalsolutions in an axisymmetric domain is given.II. D ERIVATION OF THE METHOD
The magnetostatic equations for magnetic field H ( r ) , magnetic flux density B ( r ) , and current density J ( r ) are curl H = J , (1) div B = , (2)where H is related to B via the constitutive relation H = ˆ ν B (3)with local reluctivity (inverse permeability) tensor ˆ ν ( r ) = ˆ µ ( r ) − . Written in terms of a vector potential A ( r ) with curl A = B , such that Eq. (2) is automatically fulfilled, they are equivalent to the curl-curlequation, curl ( ˆ ν curl A ) = J , (4)as long as the domain is simply connected.Our goal is to solve Eq. (4) on a finite three-dimensional domain Ω with a symmetry direction e = grad x along which the cross-section Ω t doesn’t change and where all sources and boundary conditionsare π -periodic in x . An illustration in appropriately ordered cylindrical coordinates x = Z , x = R , x = ϕ is found in Fig. 1. A. Vector calculus in 3D and 2D curvilinear coordinates
Let ( x , x , x ) be right-handed curvilinear coordinates that uniquely describe any position r = r ( x , x , x ) given by Cartesian components r k in the relevant domain. Co- and contravariant basis vectors are respec-tively defined by e k = ∂ r ∂ x k and e k = grad x k , (5)and co- and contravariant components of vector fields V ( r ) by V k = V · e k and V k = V · e k , respectively.Components g i j of the metric tensor ˆ g are g kl = e k · e l = g kl ( x , x , x ) , (6) Figure 1. Example of an axisymmetric domain with indication of the coordinate axes for x , x , x . Note that the indicated transversecross-sections do not change their shape along the symmetry direction x . and the metric determinant g = det [ g kl ] is always positive, as we are considering right-handed systems. Itssquare-root √ g is equal to the Jacobian of the transformation to Cartesian coordinates and thus the weightof the volume element. The usual symmetry condition g kl = g lk follows from the geometric definition of ˆ g in Eq. (6).For the metric-free definition of differential operators and later use in their discretised form for numericsit is useful to introduce densities of weight W to represent scalars, vectors, and tensors. A density U (denoted in calligraphic letters) of weight W for a quantity U is defined as U = √ g W U , (7)where U can represent either a scalar, or co/contravariant components of a vector or tensor. If notstated otherwise, we use the term density for the default value W = + . Usual scalars and vector/tensorcomponents correspond to the case W = . Terms such as contravariant vector density and contravariantdensity representation of a vector will be used synonymously here for easier notation, as the conceptualdifference has no practical concequence in the present context.Table I lists the choice of input and output representation for differential operators such that the Jacobian √ g is formally removed from their definition. In this way, the coordinate-independent definitions of 3Ddivergence and gradient are div U = ∂ k U k , grad U = e k ∂ k U , (8)where the divergence acts on a contravariant vector density and yields a scalar density, and the gradientacts on a scalar field and yields a covariant vector field. The 3D curl operator curl U = (cid:15) i jk e i ∂ j U k acts on a covariant vector field and yields a contravariant vector density. It contains the Levi-Civita tensor ˆ (cid:15) with contravariant density components (cid:15) i jk = for i j k being circular permutations of , − forpermutations of , and otherwise. Here and later we use the notation ∂ k = ∂ / ∂ x k and the usual convention to sum over indices appearing twice, i.e. (cid:205) k = in Eq. (8). Table IN
ATURAL INPUT AND OUTPUT FOR DIFFERENTIAL OPERATORS IN CURVILINEAR COORDINATES , DENSITIES ARE OF WEIGHT + HERE . Operator Symbol Input Output
Gradient grad scalar covariant vectorCurl curl covariant vector contravariantvector densityDivergence div contravariantvector density scalar densityTransversegradient grad t scalar covariant vectorTransversescalar curl curl t covariant vector scalar densityTransversevector curl curl t scalar contravariantvector densityTransversedivergence div t contravariantvector density scalar densityTable IIT RANSLATION BETWEEN TERMINOLOGY OF CLASSICAL TENSOR CALCULUS AND DIFFERENTIAL GEOMETRY IN DIMENSION N . Tensor calculus Differential geometry scalar -form / scalarscalar density N -formcontravariant vector vectorcovariant vector -form / covectorcontravariant vector density ( N − ) -formrank-2 tensor densities Hodge operators We want to treat the symmetry direction x separately, splitting the problem into a longitudinal partalong the symmetry direction, and a transverse part containing the remaining two dimensions. For thenotation of two-dimensional vectors we use lowercase bold letters, i.e. u = U e + U e = U e + U e . We define two-dimensional transverse divergence and gradient operators analogous to the 3D case with div t u = ∂ U + ∂ U , grad t U = e ∂ U + e ∂ U . (9)For the curl, in contrast to 3D, both a vectorial and a scalar transverse curl operator exist with curl t U e ∂ U − e ∂ U , u = ∂ U − ∂ U , (10)yielding vector and scalar densities, respectively. Note that all four of these operations either take a scalarcorresponding to the longitudinal part as input and give a two-dimensional vector corresponding to thetransverse part as output, or vice versa.Matching input and output of operators in Table I reflect the de Rham complex [14] describing in whichorder operators may act between different function spaces. In 3D this is H grad → H curl curl → H div div → L , (11)where application of two operators in a row yields zero, e.g. curl grad U = . This is the case for any deRham complex. The 3D de Rham diagram leads from a scalar field in H to a scalar density field in L . Since the curl operator mixes vector components and to distinguish between transverse and longitudinalparts, relation (11) breaks up into two separate ones in 2D, given by H grad t → H curl t curl t → L , (12) H curl t → H div t div t → L . (13)As in 3D, both diagrams lead from a scalar field in H to a scalar density field in L . The differencelies in the use of covariant vectors in H curl t or contravariant vector densities in H div t . These two casescan be translated into each other by rotation via the 2D Levi-Civita tensor ˆ (cid:15) t with contravariant densitycomponents given by (cid:15) kl t = (cid:18) − (cid:19) . (14)Relations between 2D differential operators can be written as curl t U = ˆ (cid:15) t grad t U , (15) curl t u = div t ( ˆ (cid:15) t u ) . (16) B. Covariant formulation of classical electrodynamics
Using differential operators in the metric-free way stated above, Maxwell’s equations act on densityrepresentations of fields listed in Table III. In SI units they are written as ∂ k D k = (cid:37), (17) (cid:15) i jk ∂ j E k = − ∂ B i ∂ t , (18) (cid:15) i jk ∂ j H k = J i + ∂ D i ∂ t , (19) ∂ k B k = , (20)in any curvilinear coordinate system. Magnetostatic equations (1-2) are reproduced in the stationary limitof (19–20) with vanishing time derivatives. The constitutive relations linking excitations D k , H k to locallinear responses E l , B l are D k = ε kl E l , (21) H k = ν kl B l , (22)with permittivity ˆ ε and reluctivity ˆ ν respectively represented by contravariant W = + density and covariant W = − density components, ε kl = √ g e k · ˆ ε e l , (23) ν kl = √ g − e k · ˆ ν e l . (24)While the field equations (17–20) remain coordinate-independent, Eqs. (21–22) contain all influencefrom the metric tensor ˆ g implicitly via the basis vectors and the Jacobian in Eqs. (23–24). In particularfor a scalar ν , components of ˆ g enter the resulting covariant density representation of ˆ ν in curvilinearcoordinates, ν kl = √ g (cid:213) i , j ∂ k r i ∂ l r j νδ i j = g kl √ g ν. (25)This means that covariant reluctivity components generated by a scalar ν inherit their symmetry propertiesfrom g kl . If physical components of the permeability tensor are already given in the desired coordinate Table IIIC
ONVENTIONS TO REPRESENT ELECTROMAGNETIC SCALAR , VECTOR AND TENSOR FIELDS BY DENSITIES OF VARYING WEIGHT . Quantity Symbol Representation Weight
Metric tensor g kl covariant Inverse metric g kl contravariant Jacobian √ g scalar + Levi-Civita tensor (cid:15) ijk contravariant + Charge density (cid:37) scalar + Current density J k contravariant + Electric field E k covariant Magnetic flux density B k contravariant + Electric displacement D k contravariant + Magnetic field H k covariant Scalar potential Φ scalar Vector potential A k covariant Permittivity ε kl contravariant + Permeabiliy µ kl contravariant + Reluctivity ν kl covariant − Longitudinal reluctivity ν scalar − Transverse reluctivity ν t , kl covariant − Mod. transv. reluctivity ¯ ν kl t contravariant + Transverse Levi-Civita tensor (cid:15) kl t contravariant + frame as a matrix [ µ ( kl ) ] , covariant density components of ˆ ν can be found by taking its inverse [ ν ( kl ) ] = [ µ ( kl ) ] − and computing ν kl = √ g kk g ll √ g ν ( kl ) . (26)One can see that if ν kl are constant in certain curvilinear coordinates, physical components ν ( kl ) willusually vary locally and vice-versa. One should remark that the solution for covariant components H k incurvilinear geometry with constant physical ν ( kl ) is identical to one for Cartesian components of H witha spatially varying ν ( kl ) . Thus one could emulate curvilinear geometry for magnetostatics in flat geometryby a material with locally varying permeability, and vice versa. C. Reduction of magnetostatics to 2D by Fourier expansion
To reduce the 3D problem of Eq. (4) to a number of 2D equations we write quanties assumed to be π -periodic in x as a Fourier series f ( x , x , x ) = ∞ (cid:213) n = −∞ f n ( x , x ) e i nx . (27)The following equations concern single harmonics with n omitted as an index in the notation. To becompatible with the curl-curl equation (4) in covariant form we require a Fourier expansion of covariantvector components A k of A and contravariant vector density components J k of J .To retain linearity of terms involving x and thus avoid mode-coupling via convolution in the constitutiverelation (22), we require covariant W = − density components ν kl = ν kl ( x , x ) of the impermebilitytensor ˆ ν to be independent of the symmetry coordinate x if we allow arbitrary harmonics n for the fields. In addition, to be able to split transverse and longitudinal components of fields later, off-diagonalcomponents in x shall vanish. This means that ν kl = (cid:169)(cid:173)(cid:171) ν t , kl
000 0 ν (cid:170)(cid:174)(cid:172) , (28)with transverse reluctivity ν t , kl = (cid:18) ν ν ν ν (cid:19) (29)as a covariant W = − density and longitudinal reluctivity ν as a scalar W = − density. In the caseof a scalar permeability, those symmetry restrictions thus apply to the metric tensor and vice versa, so √ g − g kl shall be independent from x and g k = g k = for k (cid:44) . A modified transverse reluctivity ˆ¯ ν t = − ˆ (cid:15) t ˆ ν t ˆ (cid:15) t (30)is useful to introduce, with contravariant density components ¯ ν kl t = (cid:18) ν − ν − ν ν (cid:19) , proportional to the transpose of the inverse of (29). For a scalar ν and g kl of the form (28) this reducesto ¯ ν kl t = √ g g kl t g ν, (31)where g kl t is the transverse part of the (symmetric) inverse metric tensor. Under the given restrictions the non-oscillatory part n = of Eq. (4) splits into a transverse and alongitudinal part with curl t ν curl t a = j , (32) curl t ˆ ν t curl t A = J , (33)where two-component a = A e + A e are expressed as a covariant vector and j = J e + J e as acontravariant vector density. The labels transverse for Eq. (32) and longitudinal for Eq. (33) refer to A and J here. This is not to be confused with B itself, for which exactly the opposite is the case, i.e. thelongitudinal B follows from Eq. (32) and the transverse b from Eq. (33), corresponding to the first stepin Eq. (12) and (13), respectively. Using the relations between 2D differential operators in Eqs. (15-16)together with the modified transverse reluctivity ˆ¯ ν t of Eq. (30), we can re-write Eq. (33) as − div t ( ˆ¯ ν t grad t A ) = J . (34)For the oscillatory part with n (cid:44) , we set the covariant component A to zero by a gauge transformation A → A − grad ∫ x x A d x (cid:48) . (35)This corresponds to setting A → A − grad A /( i n ) for single harmonic A = A n ( x , x ) e i nx with harmonicindex n . Independent gauging in that manner is possible for each n (cid:44) . Due to the superposition principle,we can use a fully transverse vector potential with a different gauge for each individual harmonic and If instead the inverse of the transverse part of ˆ g were used, the result would be identical, e.g. for cylindrical coordinates ¯ ν ZZ t = ¯ ν RR t = ν / R .This inverted dependency on R compared to the usual Laplacian ∆ in cylindrical coordinates is characteristic for the Grad-Shafranov operator ∆ (cid:63) (see [11]). take a Fourier sum in the end. In that case the contravariant magnetic flux density components are givenby B = − i nA , B = i nA , B = curl t a . (36)Ampère’s law is transformed to ∂ ( ν curl t a ) + n ( ν A − ν A ) = J , (37) − ∂ ( ν curl t a ) + n (− ν A + ν A ) = J , (38) i n [ ∂ ( ν A − ν A ) + ∂ ( ν A − ν A )] = J . (39)Eqs. (37-38) contain only transverse components of A and J and can be rewritten by using the vector curl t operator, yielding the purely 2D problem for n (cid:44) as curl t ( ν curl t a ) + n ˆ¯ ν t a = j . (40)As opposed to the singular ungauged three-dimensional curl - curl equation (32), the additional termresulting from the fixed gauge makes Eq. (40) uniquely solveable, analogous to the 3D curl - curl equationwith a time-harmonic term. Like in the scalar Helmholtz equation − ∆ Φ + n Φ = ρ arising from Fourierexpansion of the Poisson equation, the positive sign of the term weighted by n leads to rapidly decayingsolutions, opposed to oscillating solutions which would typically result from the Fourier expansion of awave equation in time. Eq. (39) links longitudinal J and B components and can be written compactly as i n div t ( ˆ¯ ν t a ) = J . (41)Eq. (41) is automatically fulfilled via the divergence relation for Fourier harmonics div t j + i n J = , (42)which can be seen from applying div t to Eq. (40).Formally, by setting n = , Eq. (40) includes Eq. (32) as a special case, so we summarize the two asthe class of transverse equations (40) for arbitrary n , whereas the longitudinal equation (34) equivalentto Eq. (33) needs only to be solved for n = . D. Variational formulation in coordinate space
To find a variational formulation for numerical computations, Eqs. (34) and (40) are multiplied by atest function and integrated with weight √ g in coordinate space ( x , x ) parametrising the cross-sectionperpendicular to the symmetry direction of the original domain Ω with a 2D volume element d Ω t : = d x d x and line element on the boundary, d Γ t = (cid:113)(cid:0) d x (cid:1) + (cid:0) d x (cid:1) . (43)Volume integration of quantities F ( x , x ) over Ω and dividing the result by the range π of x results inan integral over Ω t of the density representation F = √ g F , π ∫ Ω F ( x , x ) d V = π ∫ π d x ∫ Ω t d x d x F ( x , x ) = ∫ Ω t F ( x , x ) d Ω t . (44)The variational form of the longitudinal equation (34) for the non-oscillatory part n = with scalar testfunction w ( x , x ) is ∫ Ω t ( ∂ k w ) ¯ ν kl t ( ∂ l A ) d Ω t − ∫ Γ t w n k ¯ ν kl t ( ∂ l A ) d Γ t = ∫ Ω t w J d Ω t . (45) Here n = ( n , n ) is the unit outward normal vector across the boundary line Γ t in coordinate space ( x , x ) and the implied sums are taken over k , l from to . For the special case w = a compatibility conditionfollows as − ∫ Γ t n k ¯ ν kl t ∂ l A d Γ t = ∫ J d Ω t , (46)fixing the Neumann term in Eq. (45) corresponding to the magnetic field parallel to the transverse boundaryto the total current through the surface x = const . within the domain.For the transverse equation (40), with vectorial test function w and denoting the coordinate vector alongthe boundary in counter-clockwise direction by s = ( s , s ) = − ˆ (cid:15) t n = (− n , n ) , we obtain ∫ Ω t curl t w ν curl t a d Ω t + n ∫ Ω t w k ¯ ν kl t A l d Ω t − ∫ Γ t w k s k ν curl t a d Γ t = ∫ Ω t w k J k d Ω t , (47)for the case n = as well as n (cid:44) . E. Cartesian and cylindrical coordinates
For Cartesian coordinates and scalar reluctivity ν , the variational form of Eq. (47) with Neumannboundary condition on Γ t and test function w is given by ∫ Ω t curl t w ν curl t a d x d y + n ∫ Ω t w · ν a d x d y − ∫ Γ t s · w ν curl t a d Γ = ∫ Ω t w · j d x d y . (48)Here, s = (− n y , n x ) is the tangential vector along the boundary line Γ t of the x y cut of the domain, Ω t ,and n = ( n x , n y ) the unit outward normal vector orthogonal to s in this cross-section.For cylindrical coordinates x = Z , x = R , x = ϕ , the non-vanishing metric components are g = g = and g = R and Eq. (47) with a scalar permeability ν yields ∫ Ω t ( ∂ Z w R − ∂ R w Z ) R ν ( ∂ Z A R − ∂ R A Z ) d R d Z + n ∫ Ω t ν R ( w R A R + w Z A Z ) d R d Z − ∫ Γ t ( w R s R + w Z s Z ) R ν ( ∂ Z A R − ∂ R A Z ) d Γ t = ∫ Ω t ( w R J R + w Z J Z ) d R d Z . (49)III. F INITE ELEMENT DISCRETISATION
For the general transverse variational problem in Eq. (47), the natural discretisation for a are 2D(Nédélec) edge elements conforming to H ( curl , Ω t ) . Due to its fixed divergence in Eq. (33), components J k of the transverse current density j should be discretized by 2D (Raviart-Thomas) elements conformingto H ( div , Ω t ) . For practical reasons, it is convenient to define a quantity t in H ( curl , Ω t ) with div t j = curl t t instead. Comparison of components yields T = −J , T = J , (50)which can also be expressed via the Levi-Civita tensor defined in Eq. (14) as j = ˆ (cid:15) t t . (51)The 2D vector t is related to the n th harmonic of the current vector potential T ∝ e i nx with J = curl t T ,what can be seen from the corresponding expressions. As with Eq. (35), a fully transverse T for n (cid:44) can be chosen as T = T e + T e , to which t of Eq. (14) is proportional with factor i n . For n = , onecan instead choose a longitudinal T = T e with scalar stream function T ( x , x ) and j = curl t T , as longas J = . The re-ordering from ( R , ϕ, Z ) to ( Z , R , ϕ ) is required to fit the general framework with symmetry coordinate x . IV. V
ALIDATION RESULTS
Here we present results for an implementation of the presented Fourier-FEM approach in cylindricalcoordinates. To discretize the variational formulation, 2D Raviart-Thomas (longitudinal) and Nédélec(transversal) finite elements of lowest order are used inside FreeFEM [15]. Validity and performanceof numerical compuations are assessed based on analytical test cases with field components defined ina piecewise manner over cylinder radius R . For the numerical computations only boundary values andpossible volumetric currents are imposed.For the axisymmetric ( n = ) case we introduce an analytical longitudinal magnetic field B (cid:107) e Z with B Z = R − R ( R < . ) , R ( . ≤ R < . ) , R ( R ≥ . ) , (52)resulting in J R = , J Z = , (53)and J ϕ = J ( ϕ ) = (cid:40) ( R < . ) , ( R ≥ . ) , (54)and a second case with a transverse field B (cid:107) e ϕ given by B ϕ = R ( R < . ) , R ( . ≤ R < . ) , R ( R ≥ . ) , (55)resulting in J R = , J ϕ = , (56)and J Z = RJ ( Z ) = (cid:40) R ( R < . ) , ( R ≥ . ) . (57)As an analytical test case for the non-axisymmetric Fourier harmonic n = , we use a radial transversefield with J = and B Rn = R ( R < . ) , R − R ( . ≤ R < . ) , R + R ( R ≥ . ) . (58)Fig. 2 shows the convergence of numerical computations. The convergence rate of the L2 error overdegrees of freedom is linear, as expected from the discretization in the lowest order edge element space.V. C ONCLUSION AND O UTLOOK
An approach for the numerical solution for the magnetostatic field on 3D domains with a symmetrydirection has been described. In particular it is applicable to translationally symmetric and axisymmetricdomains. Its validity has been demonstrated on model problems in cylindrical coordinatesSince eachharmonic is computed separately, batch parallel computations are easily possible using the describedmethod.Even though arbitrary spatial variations of current density and boundary conditions can be treated, theapproach is most efficient for symmetric current distributions reducing the number of non-zero harmonics.In axisymmetric domains this includes circular ( n = ) or square-like shapes ( n = , , , . . . ). Furthermore,smoother dependencies on the symmetry coordinate lead to a faster decay in the spectrum. A more severe degrees of freedom − − − − r e l a ti v ee rr o r( L ) n=0 longitudinaln=0 transversen=1 transverse Figure 2. Convergence of numerical Fourier-FEM computations for analytical test cases. As expected the relative error decreases linearlywith the number of degrees of freedom. limitation for engineering applications, apparent in the chosen model problems, is the restriction on theshape of regions with different permeability following the symmetry direction. To lift this requirement,coupling of different harmonics would have to be taken into account, thereby removing the originaladvantage of trivial parallelisability.The method may be generalised further by considering an expansion in a different set of functions, e.g.Bessel functions for the expansion in cylindrical harmonics in the radial direction. For the treatment oftime-dependent and time-harmonic eddy current or full electromagnetic wave problems the introductionof a an additional electromagnetic scalar potential will be required for the oscillatory harmonics if thegauge is fixed to eliminate one component of the vector potential.A
CKNOWLEDGEMENTS
The authors would like to thank Fabian Weissenbacher for supporting initial investigations. Specialthanks go to Friedrich Hehl for insightful discussions on the covariant formulation of electromagnetism.The authors gratefully acknowledge support from NAWI Graz, from the OeAD under the WTZ grantagreement with Ukraine No UA 04/2017 and from the Reduced Complexity Models grant number ZT-I-0010 by the Helmholtz Association of German Research Centers.R
EFERENCES [1] J. G. V. Bladel,
Electromagnetic Fields . John Wiley & Sons, Inc., 2007.[2] J.-M. Jin,
The finite element method in electromagnetics . John Wiley & Sons, 2015.[3] C. Bernardi, M. Dauge, Y. Maday, and M. Azaiez,
Spectral methods for axisymmetric domains .Gauthier-Villars Paris, 1999, vol. 3.[4] P. Lacoste, “Solution of Maxwell equation in axisymmetric geometry by Fourier series decompostionand by use of H(rot) conforming finite element,”
Numer. Math. , vol. 84, no. 4, pp. 577–609, 2000.[5] M. Oh, “de Rham complexes arising from Fourier finite element methods in axisymmetric domains,”
Computers & Mathematics with Applications , vol. 70, no. 8, pp. 2063–2073, 2015.[6] M. F. Heyn, I. B. Ivanov, S. V. Kasilov, W. Kernbichler, I. Joseph, R. A. Moyer, and A. M. Runov,“Kinetic estimate of the shielding of resonant magnetic field perturbations by the plasma in DIII-D,”
Nuclear Fusion , vol. 48, no. 2, p. 024005, 2008.[7] B. M. Dillon and J. P. Webb, “A comparison of formulations for the vector finite element analysisof waveguides,”
IEEE Trans. Microw. Theory Tech. , vol. 42, no. 2, pp. 308–316, 1994.[8] J. A. Schouten,
Tensor analysis for physicists , 2nd ed. Dover, 1990.[9] E. J. Post,
Formal structure of electromagnetics: general covariance and electromagnetics . Dover,1997. [10] F. Gronwald, F. W. Hehl, and J. Nitsch, “Axiomatics of classical electrodynamics and its relationto gauge field theory,” Physics Notes, assembled by Dr. Carl E. Baum. , no. 14, 2005. [Online].Available: https://arxiv.org/pdf/physics/0506219[11] W. D. D’haeseleer, W. N. G. Hitchon, J. D. Callen, and J. L. Shohet,
Flux coordinates and magneticfield structure . Springer Berlin Heidelberg, 1991.[12] T. Garavaglia and J. Gomatam, “The Schrödinger equation in helical coordinates,”
Annals of Physics ,vol. 89, no. 1, p. 1, 1975.[13] O. Bíró, “Edge element formulations of eddy current problems,”
Computer Methods in AppliedMechanics and Engineering , vol. 169, no. 3-4, p. 391, 1999.[14] F. Brezzi and M. Fortin,
Mixed and hybrid finite element methods . Springer Science & BusinessMedia, 2012, vol. 15.[15] F. Hecht, “New development in FreeFem++,”