ABSTRACT (original article): This study provides analytic solutions for the magnetic field of coils and magnets that have a non-axisymmetric cylindrical geometry with a rectangular cross-section. New analytic solutions are provided for radially magnetized permanent magnet arcs, thin coil disc sectors, and thick coil sectors. If components of the 3D field are not representable in closed-form or as canonical Legendre elliptic integrals, the exact solution is given in terms of a series of regularized beta functions. The limit and hence spatial convergence is found to these series, giving a well-defined and fast solving algorithm for computation. The equations can be readily applied to find the magnetostatic field in linear or non-linear systems that contain a large set of elements. Example applications are provided to demonstrate how the field can be used to calculate forces and benchmark computational efficiency of the equations. A thorough review of the preceding literature and background theory is provided before a detailed methodology obtaining the analytic solutions contained in this compendium, and further related geometries in cylindrical or spherical coordinates. This is the first study to comprehensively solve the field equations for this collection of electromagnetic geometries. CITATION (original article): Forbes, M., Robertson, W.S.P., Zander, A.C. and Paulides, J.J.H. (2024), The Magnetic Field from Cylindrical Arc Coils and Magnets: A Compendium with New Analytic Solutions for Radial Magnetization and Azimuthal Current. Adv. Phys. Res., 3: 2300136. https://doi.org/10.1002/apxr.202300136

Introduction

This Mathematica notebook is supplementary material for the research article ‘The Magnetic Field from Cylindrical Arc Coils and Magnets: A Compendium with New Analytic Solutions for Radial Magnetisation and Azimuthal Current’ ( https://doi.org/10.1002/apxr.202300136 ). The notebook contains all analytic solutions from the article in an electronic form. The purpose is to aid reproducibility of equations (in the unfortunate event of a typographical error) and provide validation test cases for the reader. The reader should have confidence in the results before adopting them. The notebook requires a Nomenclature package, that is consistent with the research article notation (found in the Appendix).
Basic wrapper functions for the analytic solutions are provided to tabulate results at specific field points, along with a sanity check comparison to the original numeric integral. All tabulated results were compared to an FEA model within the research article. Reviewers of the research article were interested in a basic computational benchmark of the analytic solution results compared to the numerical integral solution and FEA. Partial code to this effect is included. Further to the tables, some example plots are provided to visualise: the geometry; the location of the chosen field points; the magnitude and direction of the field passing through the geometry.

0.1 Notes

◼
  • A related online repository can be found at: https://github.com/AUMAG/mag-cyl-field.
    • ◼
  • The code below is not written to optimise computational speed, but to be readable, portable, and correct.
    • ◼
  • For simplicity, the wrapper functions provided do not handle list inputs for ρ, φ, z.
    • ◼
  • Similarly, algorithms for computing a partial sum of β,δ (η,ζ,ι) are omitted and functions take a direct input for the number of terms. The number of terms are identical between all series at one field point; therefore, we only evaluate the volume coil (Section 2.3) for P=150 terms to reduce computation time. Mathematica is holding these terms symbolically before converting to a numeric value, where there ends up being a large number of ‘0.’ terms truncated to $MachinePrecision. This is not an issue with an algorithm and adaptive convergence parameters.
    • ◼
  • Mathematica uses complex numbers for the evaluation of ArcTan, ArcTanh,... and as such these return real numbers of the form a+0. when converted to a numerical value using '//N'. As such we additionally use '//Chop' to remove numerical 0. terms.
    • ◼
  • We also use ‘//Chop’ for the comparison of the analytic and numeric result, as dependent on the numeric procedure, we only expect accuracy to a certain precision.
    • ◼
  • Function handles are called from ‘ResultTable[]’ to return either the cylindrical or Cartesian (on axis) field components at a particular point.
    • ◼
  • 0.3 Geometry, 0.4 Constants, and 0.5 Field Points can be freely changed. Default values conform with those of the article cited in 0.1.
    • ◼
  • In this first revision, section 0.8 has been included to show some basic example plots using the field data, with updates to 0.3, 0.4, and 0.7.

    • 0.2 Import Variables and Functions

    SetDirectory[NotebookDirectory[]];​​<< Nomenclature` ​​<< Carlson` (*Used for some backwards compatibility with a special case of EllipticF[ϕ,1], and not required in 12.3 or later: https://reference.wolfram.com/language/ref/CarlsonRC.html https://github.com/tpfto/Carlson/tree/v1.1 *)​​IntOptions ={WorkingPrecision$MachinePrecision, AccuracyGoal53Log[10,2], PrecisionGoal24Log[10,2]};

    0.3 Test Geometry

    In[]:=
    ρʹ = 
    3
    1000
    ,
    8
    1000
    ; φʹ = -
    π
    6
    ,
    3π
    5
    ; zʹ = 
    1
    1000
    ,
    5
    1000
    ; (*Definite integral limits*)
    Visualisation of the chosen geometry.
    HollowCylinderSector[{ρ1_,ρ2_},{φ1_,φ2_},{z1_,z2_},PlotOptions_]:={​​Table[ParametricPlot3D[{r Cos[t],r Sin[t],z},{z,z1,z2},{t,φ1,φ2},Evaluate[PlotOptions]],{r,{ρ1,ρ2}}],​​Table[ParametricPlot3D[{r Cos[t],r Sin[t],z},{r,ρ1,ρ2},{t,φ1,φ2},Evaluate[PlotOptions]],{z,{z1,z2}}],​​Table[ParametricPlot3D[{r Cos[t],r Sin[t],z},{r,ρ1,ρ2},{z,z1,z2},Evaluate[PlotOptions]],{t,{φ1,φ2}}]}​​PlotOptions ={MeshNone,BoundaryStyleOff,PlotStyle{Gray,Opacity[1]}};​​Show[HollowCylinderSector[ρʹ,φʹ,zʹ,PlotOptions],BoxedFalse,AxesFalse]
    Out[]=

    0.4 Constants

    In[]:=
    u0 = 4π*
    -7
    10
    ;​​M = 955*
    3
    10
    ; (*A/m*)​​Iφ = 20; (*A*)​​Kφ = 4*
    4
    10
    ; (*A/m*)​​Jφ = 1*
    6
    10
    ; (*A/
    2
    m
    *)​​φ☆ =
    π
    6
    ; (*Diametric magnetisation direction*)​​P = 150;(*Number of terms in partial sum*)

    0.4 Test Field Points

    Seven field points to test all equations in article.
    In[]:=
    ρ1 =
    9
    1000
    ; φ1 =
    5π
    24
    ; z1 =
    31
    10000
    ; (*Standard*) ​​ρ2 =
    7
    1000
    ; φ2 = φ1; z2 = z1; (*Inside magnet/coil*) ​​ρ3 = 0; φ3 = 0; z3 = zʹ[[2]]; (*On magnet/coil axis*)​​ρ4 = ρ1; φ4 = φʹ[[1]]; z4 = zʹ[[2]]; (*Singular plane 1*)​​ρ5 = ρʹ[[2]]; φ5 = -φ1; z5 = zʹ[[2]]; (*Singular plane 2*)​​ρ6 = ρʹ[[2]]; φ6 = φʹ[[1]]; z6 =
    6
    1000
    ; (*Singular plane 3*)​​ρ7 =
    2
    1000
    ; φ7 = φ1; z7 = z1; (*Coil high field area*)
    Visualisation of the chosen field points.
    In[]:=
    PlotOptions ={MeshNone,BoundaryStyleOff,PlotStyle{Black,Opacity[.05]}};​​FieldPoints = {{#[[1]] Cos[#[[2]]],#[[1]] Sin[#[[2]]],#[[3]]}}&/@{{ρ1,φ1,z1},{ρ2,φ2,z2},{ρ3,φ3,z3},{ρ4,φ4,z4},{ρ5,φ5,z5},{ρ6,φ6,z6},{ρ7,φ7,z7}};​​ShowListPointPlot3D[FieldPoints,BoxedFalse,AxesFalse,PlotStyleColorData["Rainbow"]/@Subdivide[6]],​​ HollowCylinderSector[ρʹ,φʹ,zʹ,PlotOptions],​​ BoxedTrue,AxesTrue,PlotRange{{-ρ1,ρ1},{-ρ1,ρ1},{0,z6}},FaceGrids All,AxesLabel{"x","y","z"}​​ 
    Out[]=

    0.5 Tabular Results Format

    MagCylFieldMd_,Mρ_,Mφ_,Mz_,φ☆_,P_,ρʹ_,ρ_,φʹ_,φ_,zʹ_,z_ := ​​Module{},​​IfMd ≠ 0,(*Diametric*)​​ResultTableMd,φ☆,P,ρʹ,ρ,φʹ,φ,zʹ,z,"M⟂",BdAna,BdNum,BdAnaAxis,BdNumAxis;​​;​​IfMρ ≠ 0,(*Radial*)​​ResultTableMρ,φ☆,P,ρʹ,ρ,φʹ,φ,zʹ,z,"Mρ",BρAna,BρNum,BρAnaAxis,BρNumAxis;​​;​​IfMφ ≠ 0,(*Azimuthal*)​​ResultTableMφ,φ☆,P,ρʹ,ρ,φʹ,φ,zʹ,z,"Mφ",BφAna,BφNum,BφAnaAxis,BφNumAxis;​​; ​​IfMz ≠ 0,(*Azimuthal*)​​ResultTableMz,φ☆,P,ρʹ,ρ,φʹ,φ,zʹ,z,"Mz",BzAna,BzNum,BzAnaAxis,BzNumAxis;​​;​​Null​​​​ ​​CoilFilamentField[I_,ρʹ_,ρ_,φʹ_,φ_,zʹ_,z_]:= ResultTable[I,0,0,ρʹ[[2]],ρ,φʹ,φ,zʹ[[1]],z,"Iφ",BiAna,BiNum,BiAnaAxis,BiNumAxis];​​CoilDiscField[K_,P_,ρʹ_,ρ_,φʹ_,φ_,zʹ_,z_] := ResultTable[K,0,P,ρʹ,ρ,φʹ,φ,zʹ[[1]],z,"Kφ",BkAna,BkNum,BkAnaAxis,BkNumAxis];​​CoilShellField[K_,ρʹ_,ρ_,φʹ_,φ_,zʹ_,z_] := ResultTable[K,0,0,ρʹ[[2]],ρ,φʹ,φ,zʹ,z,"Kφ",BsAna,BsNum,BsAnaAxis,BsNumAxis];​​CoilCylField [J_,P_,ρʹ_,ρ_,φʹ_,φ_,zʹ_,z_] := ResultTable[J,0,P,ρʹ,ρ,φʹ,φ,zʹ,z,"Jφ",BcAna,BcNum,BcAnaAxis,BcNumAxis];​​​​ResultTableM_,φ☆_,P_,ρʹ_,ρ_,φʹ_,φ_,zʹ_,z_,mag_,anaFn_,numFn_,anaFnAxis_,numFnAxis_ := ​​ Module{Bana,Bnum,Bcom,heading,tab},​​ Ifρ0,​​ Bana = NanaFnAxisM,φ☆,ρʹ,φʹ,zʹ,z,$MachinePrecision //Chop;​​ Bnum = numFnAxisM,φ☆,ρʹ,φʹ,zʹ,z;​​ Bcom = Chop[Bana-Bnum,
    -8
    10
    ];​​ heading = {"Bx","By","Bz"} (*Axis, Cartesian*), ​​ Bana = NanaFnM,φ☆,P,ρʹ,ρ,φʹ,φ,zʹ,z,$MachinePrecision //Chop;​​ Bnum = numFnM,φ☆,ρʹ,ρ,φʹ,φ,zʹ,z;​​ Bcom = Chop[Bana-Bnum,
    -8
    10
    ];​​ heading = {"Bρ","Bφ","Bz"} (*Cylindrical*); ​​ ; ​​ tab = TableForm[{Bana,Bnum,Bcom}, TableHeadings -> {{"Analytic", "Numeric", "Comparison 8dp"}, heading}]; ​​ CellPrint@ExpressionCell[RawBoxes[ToBoxes[tab] /. (x : GridBox[___]) :> MapAt[mag &, x, {1, 1, 1}]], "Output"]​​ ​​ ​​ResultTableForce[M_,Mp_,ρʹ_,ρ_,φʹ_,φ_,zʹ_,z_,V_,P_,U_,O_,anaFn_,numFn_] := ​​ Module{Fana,Fnum,Fcom,heading,tab},​​ Fana = N[anaFn[M,Mp,ρʹ,ρ,φʹ,φ,zʹ,z,V,P,U,O],$MachinePrecision] //Chop;​​ Fnum = numFn[M,Mp,ρʹ,ρ,φʹ,φ,zʹ,z];​​ Fcom = Chop[Fana-Fnum,
    -8
    10
    ];​​ heading = {"Fx","Fy","Fz"} ; ​​ tab = TableForm[{Fana,Fnum,Fcom}, TableHeadings -> {{"Analytic", "Numeric", "Comparison 8dp"}, heading}]; ​​ CellPrint@ExpressionCell[tab, "Output"]​​ 

    0.6 Example field solutions

    Evaluate the same magnet geometry, at the same field point, with equal magnetisation in the diametric, radial, azimuthal, and axial directions. Creates a table for each magnetisation direction and compares to 8 decimal places the analytic solution and numeric integral.

    0.7 Timed solutions

    Figure 20 from the article, generated in-part using the functions below. Requires enabling dynamics.
    Altered tables to those in Section 0.5, including RepeatedTiming[] and an additional table to show computational efficiency.

    0.8 Example Plots

    Figure 1 from the article, generated in-part using the functions below.
    Surface geometry, mesh, and visualisation.
    Calculate magnetic field.
    Organise data and regions for plotting.
    Cylindrical contour plot and texture.
    Cylindrical streamlines.
    Cartesian streamlines.
    Outputs from the above functions.

    1.0 Permanent Magnets

    1.1 Diametric Magnetisation

    1.1.0 Equations

    Analytic and Numeric function handles. Returns B={Bρ,Bφ,Bz} or B={Bx,By,Bz} (on axis).
    Magnetisation vector for field inside magnet. Returns M={Mρ,Mφ,Mz}.
    Special cases of the geometry for the analytic function handle. Replaces BdAna[] (Cylindrical) or BdAnaAxis[] (Cartesian, on axis).
    Integrands to be solved for BdNum[] (Cylindrical) and BdNumAxis[] (Cartesian, on axis).
    Summands for BdAna[] and BdAnaAxis[].
    Singularities in the summands of BdAna[] or BdAnaAxis[].

    1.1.1 Standard - Outside Magnet

    1.1.2 Special Case a. - Inside Magnet

    2.0 Coils with Azimuthal Current Density

    Integrands to be solved for B{i,k,s,c}Num[] (Cylindrical) and B{i,k,s,c}NumAxis[] (Cartesian, on axis). Common between filament, disc, shell, volume.

    2.1 Filament

    2.1.0 Equations

    Analytic and Numeric function handles. Returns B={Bρ,Bφ,Bz} or B={Bx,By,Bz} (on axis).
    Special cases of the geometry for the analytic function handle. Replaces BiAna[] (Cylindrical) or BiAnaAxis[] (Cartesian, on axis).
    Summands for BiAna[] and BiAnaAxis[].
    Singularities in the summands of BiAna[] or BiAnaAxis[].

    2.1.1 Standard - Outside Coil Radii

    3.0 Green’s Function Integrals

    This section is simply a comparison of the integral transforms, discussed in part 3 of the article.

    3.0.0 Integrals and Analytic Solutions

    Azimuthal integral
    Radial integral
    Axial integral
    Radial surface integral
    Axial surface integral
    Volume integral

    3.0.1 Azimuthal Integral

    4.0 Magnetic field derivates

    This section has the example derivatives given in part 6.2 of the article. Integrals with respect to t are not true for all φʹ (ϕ).

    4.0.0 Derivatives and Analytic Solutions

    Elliptic integral of the first kind
    Elliptic integral of the second kind
    Elliptic integral of the third kind
    Regularised beta function

    5.0 Forces between Axially Magnetised Permanent Magnets

    5.0.0 Equations

    Analytic and Numeric function handles. Returns F={Fx,Fy,Fz}.
    Special cases of the geometry for the analytic function handle. Replaces FzAna[] or FzNum[].
    Integrands to be solved for FzNum[].
    Summands and ancillary functions for FzAna[].
    Singularities in the summands of FzAna[].

    5.0.1 An axisymmetric force

    Hollow rings.
    Solid rings (m-summations removed).

    CITE THIS NOTEBOOK

    The Magnetic Field from Cylindrical Arc Coils and Magnets: A Compendium with New Analytic Solutions for Radial Magnetisation and Azimuthal Current​
    by Matthew Forbes, Will Robertson, Anthony Zander, Johan Paulides
    Wolfram Community, STAFF PICKS, July 16, 2024
    ​https://community.wolfram.com/groups/-/m/t/3222778