Quantitative Assessment of Finite-element Models for Magnetostatic Field Calculations
QQuantitative Assessment of Finite-element Models forMagnetostatic Field Calculations
J.A. Crittenden a, ∗ a CLASSE, Cornell University, Ithaca, NY 14853, United States
Abstract
We present quantitative means for assessing the numerical accuracy of staticmagnetic field calculations in finite-element models. Our calculations use thethree-dimensional Opera simulation software suite of Dassault Syst`emes. Ourneed to assess the effects of fringe fields requires such a 3D algorithm. Whilewe do discuss and compare our approach to a method of accuracy estimationused in the Opera post-processor, our methods are generally applicable to anymodel using relaxation techniques in finite-element systems. For purposes of il-lustration, we present modeling and analysis of two types of quadrupole electro-magnets presently in operation in the south arc of the Cornell Electron StorageRing (CESR). Calculations of field multipole expansion coefficients and numer-ical deviations from Maxwell’s equations in source-free regions are discussed,with emphasis on the dependence of their accuracy on changes to the finite-element model. Successive refinement steps in the finite-element model for thenon-extraction type of CESR south arc quadrupole achieve a reduction in theRMS value of the longitudinal component of the curl vector on the magnet axisby a factor of nearly 70 from 4 . × − T/m to 6 . × − T/m, which is0.0021% of the field gradient. An accuracy of 2.4% is achieved for a dodecapolecoefficient of 2 . × − of that of the quadrupole in the longitudinal integral ofthe vertical field gradient at a radius of 1 cm. ∗ Corresponding author. Tel.: +1 6072554882
Email address: [email protected] (J.A. Crittenden)
Preprint submitted to Elsevier January 18, 2021 a r X i v : . [ phy s i c s . acc - ph ] J a n ontents1 Introduction 22 Electromagnets in the CESR South Arc 5 Accurate modeling of static magnetic fields plays an important role in thedesign, construction, commissioning and operations for particle accelerators.We have developed means for quantitative assessment of the accuracy of fieldcalculations in finite-element-based models. General treatments of the finite-element method and engineering applications can be found, for example, inRefs. [1, 2, 3]. Here we use the examples of two types of quadrupole magnetused in the CHESS-U upgrade of the Cornell Electron Storage Ring (CESR) [4].The 768-m-circumference CESR ring with its injectors is shown schematicallyin Fig. 1, including the region of the south arc upgrade. The principal in-novation of this upgrade was the replacement of the linear south interactionregion formerly occupied by the CLEO detector [5] with an arc of radius sim-ilar to that of north arcs of the ring. The elimination of the strong dipolemagnets in regions of high dispersion at the entrance and exit of the southarc region permitted a reduction in beam emittance by a factor of four, result-ing in brighter X-ray beams provided to the Cornell High-Energy SynchrotronSource end stations. The south arc region consists of six double-bend achro-mats, each of which comprises two horizontally defocusing dipole/quadrupole2 igure 1: The Cornell Electron Storage Ring and injector configuration. Six-GeV positronscirculate in the clockwise direction in the storage ring. combined-function magnets and two horizontally focusing quadrupole magnets.This structure provides for one straight section for insertion devices to provideextracted X-ray beams, as shown in Fig. 2. Figure 3 shows a photograph of halfof the double-bend achromat layout. The two quadrupoles at the entrance andexit of each achromat are designed for exact quadrupole symmetry (Type A),while the two central quadrupoles are equipped with extraction channels on theoutside of the ring (Type B). Figure 4 shows color maps of the magnetizationfield vector magnitude on the steel surfaces of the Opera models described inthis report.One motivation for this work is to assess the accuracy of the idealized ap-proximate models of our quadrupole magnets in the development and operationslattice model based on the Bmad library, [6] and to determine quantitatively itsadequacy for any given purpose. These models consist merely of values for k (or field gradient) and length, chosen such that their product is consistent withmodeling and field measurements. The assessment can be done by comparisoneither with tracking through modeled field maps (see Ref. [7, 8] for a recentexample of lattice analysis using field-map tracking) or by adding field integralmultipole terms to the Bmad lattice element. The latter method is orders ofmagnitude faster in computation speed, but its accuracy must first be verified3 igure 2: Schematic diagrams of the double-bend achromat structure. by comparison with reliable tracking through field maps. Thus knowledge ofthe numerical accuracy of multipole calculations in the model is required. Ourapproach is described in Sec. 4.1. We make use of the symmetry constraints forthe quadrupole-symmetric magnet Type A to obtain estimates for the numericaluncertainties inherent in the refined finite-element model. These can then beused to estimate the accuracy in the multipole coefficients for the quadrupolemagnet Type B. By comparing multipole expansions in the vertical field com-ponent at the center of the magnet to that of the integrated field component, wecan distinguish multipole contributions by the pole shape from those contributedby the fringe field.Our second criterion for estimating the accuracy of the models is an anal-ysis of the numerical violations of Maxwell’s equations in a volume relevant tothe particle trajectories, as described in Sec. 4.2. The divergence and the threecomponents of the curl of the electric field vector deviate from zero primarilyin the fringe field region in our models. We calculate the RMS value of each ofthese quantities along longitudinal lines as a function of horizontal position toestimate the improvement in magnitude and location of the errors as the finite-element model is refined. Since the longitudinal sums of the errors tend to besmall, we use the RMS values, since all errors contribute errors in the particletrajectory. We employ a difference-over-sum method to obtain relative devi-ations along the longitudinal coordinate, then propagate those relative errorsto obtain an uncertainty in the RMS value as a function of transverse posi-tion. This dependence of non-zero divergence and curl values on the horizontalcoordinate motivates various types of changes to the finite-element model.These diagnostic criteria are applied to our Opera models for the south arcquadrupoles in Sec. 5, showing their evolution at each step of the progressiverefinement of the finite-element algorithms used in the volumes and surfaces of4 igure 3: Photograph of half of the double-bend achromat structure, showing a combined-function magnet and two quadrupole magnets, viewed from downstream. The upstreamquadrupole magnet exhibits the extraction channel for the X-ray beam. the steel and air structures.
2. Electromagnets in the CESR South Arc
The electromagnets in the south arc were designed, built and installed byS. Barrett, A. Lyndaker and A. Temnykh in collaboration with industrial part-ners and the CLASSE technical staff.
The horizontally defocusing component of the south arc optics is provided bytwelve combined-function magnets, each spanning 74.8 mrad with an arc lengthof 2.3475 m. The vertical field component on the beam trajectory is 6.377 kG.The gap increases toward the outside of the ring, where the open geometry ofthis C-shape allows the X-ray beams from the insertion devices to exit the ring.The field gradient on the beam trajectory is 8.77 T/m. The results of a detailed3D modeling study of these magnets, including multipole analysis and trackingfor 1 GeV and 6 GeV positrons are available [9].
The CHESS-U south arc optics include 24 horizontally focusing quadrupoleelectromagnets. The dimensions and nominal magnetic field and electrical pa-rameters of these magnets are shown in Table 1.5 igure 4: Color maps of the magnetization field vector magnitudes calculated in the Operamodel on the non-extraction (Type A, left) and extraction (Type B, right) south arcquadrupole magnets at nominal field intensity. The units are Gauss. Some effects of sat-uration are evident in the ends of the pole tip.
3. The Opera Application Suite
The construction of an Opera model for an electromagnet can be categorizedin three parts: • Steel volumes. These volumes can be reduced according to the symmetryof the model. For example, in the case of the Type A quadrupole, a1/16 model is possible, and the symmetries about the XY, YZ, ZX and45-degree planes in the field map are guaranteed, since field values arecalculated only in the region of the 1/16 model.Various native geometrical shapes are provided with optimized meshing al-gorithms. Algebraic functions can be modeled as well, using the MORPHcommand. The latter can be used for quadrupole pole face shapes, as cana list of 3D points, using the WIREEDGE command. The Modeler canimport files in Spatial’s ACIS solid modeling format (.sat) files) but themeshing can be much less efficient and much more time consuming. In-deed, if the sat file is too complicated, e.g. with many very small surfaces,the mesh can fail entirely. In our present case, simplifying modificationsto an SAT file were required to achieve reliable meshing , i.e. meshingsuccess robust against the variety of modifications to the model describedin Sec. 5.The accuracy of the finite-element model and overall field calculation isvariously sensitive to the parameters of maximum cell size and cell geom- The horizontal, vertical and longitudinal coordinates are denoted X, Y and Z, respectively.The origins are at the center of the quadrupole geometry. Made available by A. Lyandaker, CLASSE We acknowledge valuable assistance from Y. Zhilichev, Dassault Syst`emes. able 1: Dimensions and nominal operating parameters used in the Opera models for theCESR south arc quadrupole magnets Parameter Type A Type Bnon-extraction extraction
Number of magnets 12 12Bore radius (cm) 2.30Steel height (cm) 47.44Steel width (cm) 47.44 56.72Steel length (cm) 37.70 33.30Length including coil (cm) 54.33 49.94Pole width (cm) 8.56Central Field Gradient (T/m) -30.36 -28.01Field Gradient Integral (T) -12.31 -10.13Good Field Region ∗ (mm) ± ∗ (%) 0 - 0.1NI per coil (Amp-turns) 6579 6033Turns per coil (4 ×
8) + 7 + 6 + 5 + 4 + 3 + 2 = 59Conductor cross section (cm x cm) 0 . × .
635 with 0.318 diameter holeConductor straight length (cm) 37.7 - 47.28 33.27 - 42.78Coil inner corner radius (cm) 1.96Conductor length per turn, avg (cm) 34.6 30.6R coil (Ω) 0.0214 0.0198L (mH) 4 x 27 4 x 24Power supply current (A) 111.5 102.3Current density (A/mm ) 2.02 1.85Voltage drop/magnet (V) 9.54 8.10Power/magnet (W) 1064 829 ∗ Defined as relative deviation from the central field gradient integral less than 0 . etry type defined for each volume. These must be chosen judiciously inthe interest of model size and computation time as well.The magnetic properties of the steel material can be provided in the formof a table of values relating magnetic flux density values to magnetic fieldvalues. The south arc quadrupoles were made of a steel similar to thelow-carbon 1010 steel. • Air volumes. The encompassing volume must be chosen large enoughthat the (automatic) outer tangential magnetic boundary conditions havenegligible effect on the field calculations in the regions of interest. Airvolumes in the regions of interest, such as those relevant to the particletrajectories, must be chosen to be of sufficient refinement for the desiredaccuracy. Since these two types of volume have large differences in cellsize and volume, buffer volumes of intermediate cell size may be requiredto aid in the meshing convergence.The comments above on maximum cell size and type in the list itemconcerning steel volumes apply to the air volumes as well. It is to be7oted, however, that the Modeler meshing algorithm intrinsically adjustsits result to the local derivatives in the steel volume. As a consequence,reducing the maximum cell size parameter does not have the third-powereffect on model size and computing time that it does in the air volumes.We will also see in Sec. 5 that adjusting cell sizes on surfaces can beeffective in improving model accuracy without impractical increases inmodel size and computation time. • Conductors. The model is required to contain the full 3D geometry ofthe conductors, even when the steel volumes take advantage of availablesymmetries. In the case of coils arranged with symmetry about the Z axis,provision is made for defining the coil geometry once and specifying thesymmetry. A wide variety of native geometries are available. In the caseof our quadrupoles, we modeled each of the four coils as seven racetracks,blocks of 8x4, 7x1, 6,1, 5x1, 4x1, 3x1, and 2x1 1/4-inch square conductorsgiving the 59 turns in trapezoidal shape of the coil.
Magnetostatic calculations are implemented in the Opera package via theTOSCA program which originated at the Rutherford Laboratory in the 1960’s.Its development and many applications are documented in the proceedings ofthe Compumag conferences, which are organized by the non-profit InternationalCompumag Society. The physical basis for the TOSCA algorithm is describedby Simkin and Trowbridge in their contribution to the 1976 Compumag pro-ceedings [10].Invocation of the default iterative TOSCA program is largely parameter-free, the memory use and computing time determined by the size of the finite-element model defined in the Modeler. The calculation has multi-core capability.There is an option to use a direct solver rather than the iterative method, withattendant greater memory requirements. The default results of the TOSCAcalculation are values for the magnetic scalar potential, magnetic field vectorand magnetic flux density vector at each node in the model. There is provisionfor running multiple field calculations with varying excitation current for thepurpose of producing saturation curves.
The Post-processor can calculate and write tables on a regular grid on thebasis of the TOSCA result in two ways. The first way is to interpolate the fieldvalue at any point by averaging over nearby nodes. The second method is tointegrate over all current density and steel magnetization sources. In our caseof calculating field values in the beam region, this latter method produces, byconstruction, small local derivatives in the result, since the sources are far awayand closer to each other compared to the table step size. Such a result can ap-pear attractive, since the residuals of a polynomial fit for multipole coefficients,for example, are smaller than those from the local node interpolation method.8owever, the overall accuracy of the field calculation depends on the mesh ac-curacy everywhere . A finer mesh at the beam improves the field calculation atthe distant nodes as well. So, perhaps counter-intuitively, refinement the localmesh at the beam can improve the accuracy of the integration method. Theintegration method, despite its smoothness, requires a sufficiently detailed meshat the beam, which obviously improves the method of nodal interpolation. Thusthe use of the integration method does not by itself improve the accuracy of thefield values in the table. Mesh tuning is necessary in either case, and since thenodal interpolation method must reach the required accuracy in any case, theintegration method can be considered superfluous. Indeed, the original moti-vation for including the integration method option was not for improved localaccuracy, but rather to provide better field calculations in volumes far from thesteel, where the mesh may be very coarse. The integration method can be misleading in another way as well. Whenboundary conditions are defined on symmetry planes in the beam region, theintegration method provides an exact constraint by definition, because there areno sources on the constrained plane. For example, the vertical field componenton the magnet axis from magnetization sources is calculated to be zero bydefinition, since the sources are only available in the 1/8 model, far from thebeam. This provides a distortion of the numerical model on the plane, withnodes neighboring the plane exhibiting larger derivatives to the plane than toother neighboring nodes. Again, the only way to improve obedience to thesymmetry condition is to refine the mesh. The method of nodal interpolationproduces violations of Maxwell’s equations on the symmetry planes consistentwith the accuracy of the finite-element model.So what criteria can be used to assess the accuracy of the field calculationsin the field map? One idea might be to use the difference between the two typesof field calculations, but the above discussion shows that such a comparison isfraught with systematics unrelated to the overall accuracy of the model. TheOpera Post-processor provides a system variable which addresses the local ac-curacy of the field calculation. It compares the nodally interpolated field valueto an average of scalar potential derivatives in the nodes making up the localcell. This variable, however, provides only a value for the numerical uncertaintyin the magnitude of the magnetic field vector. Its intended purpose is to revealregions in which the mesh requires refinement by comparison to other regions,not to provide a quantitative estimate of the field calculation accuracy.The quest for such a quantitative estimate motivated the present work. Ourapproach is described in Sec. 4 below, using the CESR south arc quadrupolemodels for illustrative purposes.——————————————————————————- Private communication from Y. Zhilichev, Dassault Syst`emes . Field Calculation Diagnostic Criteria Numerical errors in our finite-element model can result in a multipole fitgiving disallowed multiple coefficients (“phantom multipoles”) even in a per-fectly symmetric geometry, as in the case of the south arc quadrupole Type A.As discussed in Sec. 3.3 the integration method of field calculation can enforcethe symmetry at the cost of anomalous discontinuities, but for our applicationwe wish to use the violations of symmetry in our field map as a measure of thenumerical accuracy of the successive models as we refine them (see Sec. 5), so weuse the nodal interpolation method of calculating the field values. Confusion canalso arise from correlations between the coefficient values in the polynomial fit,so it is advantageous to omit poorly determined coefficients. For the purposes ofthis illustration, we omit the skew multipole coefficients and fit a polynomial tothe horizontal dependence of the longitudinal integral of the vertical field com-ponent. The field map extends over a horizontal region of -10 mm < X <
10 mmin 21 1-mm-wide steps and over a longitudinal region -300 mm < Z <
300 mm in301 2-mm-wide steps. The result of a number of trial procedures is the following:1. Perform a polynomial fit to (cid:82) B Y (X , Y = 0 , Z) dZ. The result of this stepfor the final model in our refinement procedure is shown in Fig. 5. We usethe optimization package MINUIT [11, 12] with the enhanced uncertaintyevaluation referred to as MIGRAD+HESSE+MINOS [13, 14, 15]. Thesingle value chosen for the weight on each squared residual in the χ calculation is arbitrarily set to 1 Tm − , then the fit is repeated with weightset to a value such that χ /NDF is close to unity, as shown in Fig. 5. Thecorresponding weight value in this case is σ − = (1 . × − Tm) − . Wemonitor this measure of the scatter of the field integral around the fitresult during the model refinement procedure described below in Sec. 5.2. Now consider the distribution of the residuals as shown in Fig. 6. If anysmooth periodic structure indicates the need for additional terms in thepolynomial, repeat the fit with more terms. Use the fewest terms possibleto avoid complicated correlations between the results for the coefficients. The granularity of the field map is an additional source of numerical errors, beyond thosepresented by the finite-element model. While we do not address these in this work, we notethat the diagnostic methods described here can be used to quantify their contributions aswell. Multipole Table Integrated b_n Relative to 4-pole at R=0.01 m in units of 1e-4Dipole 0 Tm 0.0002Quadrupole -12.2277 Tm/m 10000Sextupole 0 Tm/m-2 -0.0007Octupole 1.2451 Tm/m-3 -0.1018Decapole -0.7881 Tm/m-4 0.0006Dodecapole -296143 Tm/m-5 2.4219
X (m) S t r a i g h t - li n e ve r ti c a l f i e l d i n t e g r a l ( T m ) s = 1.5E-07 Tm -0.1-0.0500.050.1 -0.01 -0.0075 -0.005 -0.0025 0 0.0025 0.005 0.0075 0.01 Figure 5: Initial fit for the multipole coefficients in (cid:82) B Y (X , Y = 0 , Z) dZ. The weights in the χ calculation are chosen to give χ /NDF near unity so that the uncertainties given for thecoefficient calculations are reasonable.Figure 6: Residual distribution resulting from the fit result shown in Fig. 5. The lack of asmooth power-law structure indicates that no terms need to be added to the polynomial.
11. Subtract the result for the linear terms in the fit and fit again, allowing allcoefficients to vary. This provides a more accurate determination of thehigher multipoles. Figure 7 shows the result in the present example. Select
X (m) S t r a i g h t - li n e ve r ti c a l f i e l d i n t e g r a l ( T m ) w it h li n e a r t e r m s s u b t r a c t e d Fit result for higher multipoles s = 1.5E-07 Tm -0.3-0.2-0.100.10.20.3x 10 -4 -0.01 -0.0075 -0.005 -0.0025 0 0.0025 0.005 0.0075 0.01 Figure 7: Fit result for the higher multipoles with the linear contributions subtracted. any coefficients smaller than the uncertainty in their determination andforce them to zero. In the present example, we have an octupole term of − . × − ± . × − Tm/m and a decapole term of 0 . ± . .12. Repeat the full fit with these coefficients constrained to zero and reset theweight value to obtain χ /NDF near unity. This is a small correction.The result is shown in Fig. 8. Multipole Table Integrated b_n Relative to 4-pole at R=0.15 m in units of 1e-4Dipole 0 Tm 0Quadrupole -12.2277 Tm/m 10000Sextupole 0 Tm/m-2 0Octupole -0.151 Tm/m-3 2.7787Decapole 0 Tm/m-4 0Dodecapole -285362 Tm/m-5 118146X (m) S t r a i g h t - li n e ve r ti c a l f i e l d i n t e g r a l ( T m ) Constrained s = 1.8E-07 Tm -0.1-0.0500.050.1 -0.01 -0.0075 -0.005 -0.0025 0 0.0025 0.005 0.0075 0.01 Figure 8: Polynomial fit to the full distribution with the insignificant (or poorly determined)multipole terms removed.
13. Finally, subtract the linear terms and fit to find the definitive values for thehigher-order multipoles, as shown in Fig. 9. We find significant octupoleand dodecapole coefficients of − . ± .
89 Tm/m and 2 . ± . × Tm/m . X (m) S t r a i g h t - li n e ve r ti c a l f i e l d i n t e g r a l ( T m ) w it h li n e a r t e r m s s u b t r a c t e d Constrained fit result for higher multipolesConstrained s = 1.8E-07 Tm -0.3-0.2-0.100.10.20.3x 10 -4 -0.01 -0.0075 -0.005 -0.0025 0 0.0025 0.005 0.0075 0.01 Figure 9: Results of the fit for higher multipoles after the linear part has been subtracted andthe insignificant multipole terms have been forced to zero. B Y (X , Y = 0 , Z) The result, shownin Fig. 10, indicates dodecapole is primarily sourced by the internal pole face
X (m) V e r ti c a l f i e l d c o m po n e n t a t Y = Z = ( T ) w it h li n e a r t e r m s s u b t r a c t e d Constrained fit result for higher multipolesConstrained s = 2.87E-07 Tm -0.1-0.075-0.05-0.02500.0250.050.0750.1x 10 -3 -0.01 -0.0075 -0.005 -0.0025 0 0.0025 0.005 0.0075 0.01 Figure 10: Result of the procedure to determine the higher multipole content, now appliedto the central field B Y (X , Y = 0 , Z) rather than to its longitudinal integral. The largerdodecapole terms shows that the fringe field compensates the term in the central field, reducingthe term in the integral by a factor of nearly three. shape, with the fringe field reducing its contribution in the integral by a factorof about three. The octupole term is relatively poorly determined, preventingany conclusion. 15pplying the above procedure to the field calculation for the quadrupoleType B can now identify the effect of the extraction channel on the highermultipole content. The result for the field integral shows significant sextupoleand decapole terms in addition to the octupole and dodecapole terms, as shownin Fig. 11. The octupole and dodecapole terms are of magnitude similar to
X (m) S t r a i g h t - li n e ve r ti c a l f i e l d i n t e g r a l ( T m ) w it h li n e a r t e r m s s u b t r a c t e d Constrained fit result for higher multipolesConstrained s = 1.18E-07 Tm -0.4-0.3-0.2-0.100.1x 10 -4 -0.01 -0.0075 -0.005 -0.0025 0 0.0025 0.005 0.0075 0.01 Figure 11: Application of the procedure for determining higher multipole content to the southarc quadrupole Type B shows that the extraction channel introduces sextupole and decapoleterms in the integral of the vertical field component. those of the Type A quadrupole, though the octupole term is of opposite signand much better determined. Figure 12 shows that the sextupole term in theintegral is significant, but about a factor of four smaller than at the longitudinalcenter of the magnet. The same is true of the dodecapole term. The decapoleterm is due to the fringe field. The fringe field overcompensates the octupoleterm, flipping its sign. The small decapole term is introduced by the fringe field.16
X (m) V e r ti c a l f i e l d c o m po n e n t a t Y = Z = ( T ) w it h li n e a r t e r m s s u b t r a c t e d Constrained fit result for higher multipolesConstrained s = 3.54E-07 Tm -0.125-0.1-0.075-0.05-0.02500.0250.05x 10 -3 -0.01 -0.0075 -0.005 -0.0025 0 0.0025 0.005 0.0075 0.01 Figure 12: The result shown here for the multipole analysis to the central profile of the verticalfield component indicates that the fringe field compensates (reduces) both the sextupole anddodecapole terms caused by the extraction channel at the center of the magnet by about afactor of four.
17o summarize this section on multipole analysis, we have used the exam-ple of the final, fully refined, model to exemplify the accuracy with which themultipoles in the field integral have been determined. Multipole coefficients aredetermined at the level of a few percent. The field value fluctuations aroundthe polynomial fit are about 10 − Tm, which is about 10 − of the field integralat x = 0 .
01 m. Section 5 will describe how the various steps in the refinementof the finite element model lead to this result. We note that the required degreeof refinement depends crucially on the magnet geometry. In our present case ofa well-defined pole tip shape and a magnet length much greater than the boreradius, we require a highly accurate field calculation to determine the smallhigher multipole contributions. 18 .2. Numerical Violations of Maxwell’s Equations
The second means of quantifying numerical sources of error in the field cal-culation to be illustrated is the deviation from Maxwell’s equations. Since thereare no sources in the region of interest to particle tracking, the values of thedivergence and curl should be small in an accurate finite-element model. Thevalues for the derivatives of the field components are calculated for each pointin the field map using the two adjacent points in each of the three coordinates.Edge values are omitted in the following plots. Again using the example of thefinal refined model, we show in Fig. 13 the values of the divergence and the curl D i ve r g e n ce ( T / m ) o n ax i s C u r l X ( T / m ) o n ax i s C u r l Y ( T / m ) o n ax i s C u r l Z o n ax i s ( T / m ) Z (m) a)b)c)d)
Figure 13: Values of a) the divergence and the b) X, c) Y, and d) Z curl components asfunctions of Z on the magnet axis for the quadrupole Type A. While the divergence is accurateat machine accuracy, the curl components exhibit errors in the fringe field regions. components as a function of the longitudinal coordinate Z on the magnet axisfor the quadrupole Type A. The divergence values are accurate at a level simi-lar to machine accuracy, but the curl values show errors at the 5 × − T/mlevel in the fringe field regions. For any given application, any of these fourerror sources may be of interest, but for our present purposes of illustration, weconcentrate on the Z component of the curl. Figure 14 shows the longitudinaldependence of the difference, sum and ratio of the vertical and horizontal gra-dients. Difference is the Z component of the curl, and the sum is close to twicethe field gradient. The relative errors reach 0.03% in the exterior of the magnet.It should be noted that this is not a general characteristic of all models, or evenof this one off axis. For example, the relative errors in the Z component of the19 C u r l Z | o n ax i s ( T / m ) | d B Y / d X + B X / d Y | o n ax i s ( T / m ) Z (m) R a ti o o n ax i s ( % ) dB x /dy - dB y /dxdB x /dy + dB y /dxdB x /dy - dB y /dxdB x /dy + dB y /dxa)b)c) Figure 14: Absolute values of the a) difference, b) sum and c) ratio of the vertical andhorizontal field derivatives as functions of the longitudinal coordinate Z. curl at X = 1 cm in the body of the magnet are comparable to those outside ofthe magnet, reaching a level of about 1%.20aving obtained a distribution in Z of relative errors for all X positions inthe field map, we can plot this measure of field error as a function of X. Wechoose to use the longitudinal sum of squared errors, since each contributes anerror in particle tracking, regardless of sign. For this reason, the simple sumwould provide a misleading measure due to cancellations. Also, an uncertaintyin the sum of squares of relative errors is readily obtained, providing a valuefor the uncertainty in the error obtained for each X value. The results for thesummed squares and their square root as functions of X are shown in Fig. 15.The vertical scale in Fig. 15 a) is in units of 10 − to aid the comparison to Su mm e d s q u a r e s o f r e l a ti ve e rr o rs o n Y = p l a n e ( × ) R oo t o f s u mm e d s q u a r e s o n Y = p l a n e ( % ) X (m) a)b) -4 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.00800.050.10.150.20.250.30.350.4x 10 -2 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 Figure 15: Results for a) the sum over Z of the squared ratios of the difference of dB Y /d Xand dB X /d Y to their sum, and b) the square root of the summed squares, as a function ofthe horizontal coordinate X.
Fig. 15 b), where the vertical scale is in percent.Thus we have obtained a quantitative measure of the field calculation errorswhich can be compared from model to model during the refinement of the finite-element mesh. Such a comparison is shown in Fig. 16, where the refined version 6is compared to the initial version 1 (see Sec. 5 below for elaboration.) One canconclude that the refinement process emphasized the region of interest near thebeam on the magnet axis and reduced the relative error by a factor of more Tracking results suffer in addition the accumulation of errors along the trajectory, whichwe do not account for here. oo t o f s u mm e d s q u a r e s o n Y = p l a n e ( % ) a) R oo t o f s u mm e d s q u a r e s o n Y = p l a n e ( % ) X (m) b) Figure 16: The square root of the summed squared ratios of the difference of dB Y /d X and dB X /d Y to its sum as a function of the horizontal coordinate X for a) the initial model, andb) the final finite-element model version 6. See Sec. 5 for details on the successive refinementsof the finite-element model. One can conclude that the field errors near the magnet axis havebeen reduced by a factor of more than 70. than 70.In order to compare the progress in multipole finite-element models of pro-gressive refinement, we wish to show these results in a single figure, as shown inFig. 17 a). However, this depiction requires a vertical scale which emphasizesthe worst model, while we prefer the figure to clearly characterize the featuresof the improved model. For this reason, we choose to plot the inverse of theroot of the summed squares of the relative errors, as shown in Fig. 17 b), in ourdiscussion of model development below in Sec. 5.22 ersion 1 Quad Type A Map type nodalVersion 6 Quad Type A Map type nodal R oo t o f s u mm e d s q u a r e s ( % ) a) Version 1 Quad Type A Map type nodalVersion 6 Quad Type A Map type nodal
X (m) I n ve rs e o f r oo t o f s u mm e d s q u a r e s ( × - ) b) Figure 17: a) Comparison of the results for the root of the summed squares of field errors for thetwo models shown in Fig. 16. This comparison has the distinct disadvantage of emphasizingthe worse model over the better one. In order to clearly show the features of the more refinedmodels, we will use the inverse of these values in our discussion of model development to followbelow, as shown in b). . Finite-element Model Development We begin our review of the development of the finite-element model for theCESR south arc quadrupoles with a description of the initial model, which wecall Version 1. This model takes advantage of the full quadrupole symmetryof type A, so it is a 1/16 model as shown in Fig. 18. This geometry allows
Figure 18: Geometry and finite-element mesh of the initial model (version 1). Note the1/16 volumes for the air and steel, taking full advantage of the quadrupole symmetry and thesymmetries about the XY, XZ, YZ and 45-degree planes. for tangential magnetic boundary conditions on the XY and 45-degree surfacesof the steel and air volumes, and for normal magnetic boundary conditions onthe XZ and YZ surfaces. Three air volumes are defined: 1) a small volumearound the magnet axis, 2) a buffer volume providing transition between thefine mesh of the beam volume and the other volumes, 3) a background volumewith a coarse mesh (not shown in Fig. 18) which defines the distance from theinner volumes to the tangential boundary conditions on its outer surfaces. Thebuffer volume is defined as a cylinder with the pole steel cut away. There aretwo steel volumes: 1) the intersection of a cylinder with the pole tip to allowthe definition of a finer mesh in the pole tip than in the rest of the steel, and2) the rest of the magnet steel. The dimensions, maximum cell size specificationand cell geometry type for each volume of the full model are listed in Table 5.1.Maximum cell size parameters can be defined both in volumes and on surfaces.24 a b l e : F i n i t e - e l e m e n t m e s hp a r a m e t e r s o f t h e fi v e v o l u m e s i n t h e i n i t i a l m o d e l f o r t h e s y mm e t r i c T y p e A q u a d r up o l e . A ll d i m e n s i o n s a r e g i v e n f o r t h e f u ll m o d e l. V o l u m e W i d t h H e i g h t L e n g t h R a d i u s C e ll g e o m e t r y t y p e M a x i m u m ce ll s i ze ( m )( m )( m )( m )( m ) V o l u m e XY p l a n e X Z p l a n e B e a m . . . - Q u a d r a t i c . B u ff e r -- . . L i n e a r . B a c k g r o und . . . - L i n e a r . P o l e s t ee l -- . - Q u a d r a t i c . . - Q u a d r up o l e y o k e . . . - L i n e a r . .2. Refinement Steps Intermediate results for five refinement steps of the finite-element model wereobtained in order to understand the effectiveness and performance cost of eachtype of refinement.
Version 2
The maximum mesh size on the XZ and 45-degree surfaces of the beamvolume was reduced from 2.0 mm to 0.5 mm. Table 5.3.1 below in Sec. 5.3shows that the model increased in size by an order of magnitude. Decreasingthe maximum mesh size in the full beam volume would have been impractical.The accuracy of the field calculation improved remarkably, reducing the measureof deviations from Maxwell’s equations by a factor of 20, as shown in Fig. 19.However, the accuracy in the determination of the multipole coefficients is stillnot stable, as shown by their fluctuations during further refinements of thefinite-element mesh.
Version 3
Again modifying maximum cell sizes on surfaces only, this time the maximumcell size on the pole tip surfaces was reduced to 2 mm from 5 mm. Quadraticfinite-element cell shape definitions were also introduced on these surfaces. Thesmall changes in the diagnostic criteria show that the mesh was already suffi-cient.
Version 4
The maximum cell size in the beam volume was reduced from 2 mm to1 mm. The maximum cell size in the buffer volume was reduced to 5 mm from10 mm and the cell type was changed from linear to quadratic. The dramaticreduction in the fluctuations of the field values around the polynomial fit forthe multipoles ( σ Residuals ) is reduced by more than a factor of 20) shows thatrefinement of the mesh in the volume between the beam volume and the poletips is important for the overall accuracy of the model.
Version 5
An additional micro-beam surface on the XZ plane was added, just 4 mmwide and 600 mm long, with a cell size of only 0.2 mm. This reduced thecurl values on the beam axis by about a factor of three. Interpretation of theresults for the multipoles and σ Residuals are complicated by a fluctuation in thefit which caused the octupole term to be eliminated. However, the stability inthe determination of the dodecapole term is encouraging. It is determined withan accuracy of about 3%. 26 ersion 6
The final refinement step was a test of the accuracy of the fringe field con-tribution. A 100-mm-long, 50-mm wide XZ surface with a cell size of 0.4 mmlongitudinally centered on the end of the magnet steel was introduced. Thestability of the values for the integrated multipoles within their uncertaintiesindicates that their determination is not limited by the field calculation accuracyin the fringe field region.
The results of the multipole analysis described in Sec. 4.1 for each of the sixfinite-element models are shown in the Table 5.3.1, together with the CPU timerequired for the TOSCA calculation and the number of finite-element nodesin the model. A measure of the scatter of the field calculation around thepolynomial fit result, σ Residuals is taken from the full fit following eliminationof the insignificant multipole terms (see Fig. 8 for the case of the fully refinedmodel, version 6). The final results for the multipoles are taken from the laststep of the procedure, the fit to the field values with the linear terms removed andthe disallowed multipole terms constrained to zero. Since we are interested inthe sensitivity to disallowed multipoles, we also include the uncertainty in theirdetermination from the initial fit with linear terms removed prior to settingthem to zero (as shown in Fig. 5 for version 6). This is the fit result usedto determine which terms to constrain to zero. The accuracy obtained in therefinement process for the limit on the sextupole term b (decapole term b ) is3 . × − Tm/m (3 . × Tm/m ). In the convention sometimes used forthe expression of multipole coefficients relative to the main quadrupole term,these values both correspond to 0.03 units. The octupole (dodecapole) termsare found to be − . ± .
89 Tm/m (2 . ± . × Tm/m ), or 1 . ± . × − (2 . ± . We acknowledge that magnet fabrication errors at this level may be impractical and thattheir contribution to lattice errors may thus be larger. Our goal is the quantitative assessmentof errors in finite-element models. We present a method to ensure that design model errorsdo not exceed the fabrication errors. a b l e : R e s u l t s o f t h e m u l t i p o l e a n a l y s i s p r o ce du r e a pp li e d t o t h e fi e l d i n t e g r a l s (cid:82) B Y ( X , Y = , Z ) d Z f o r t h e s i x fin i t e - e l e m e n t m o d e l s . T h e v a l u e f o r σ R e s i du a l s i s t a k e n f r o m t h e f u ll p o l y n o m i a l fi t a f t e r t h e i n s i g n i fi c a n tt e r m s h a v e b ee n c o n s t r a i n e d t o ze r o . T h e un ce r t a i n t i e s i n t h e s e r e m o v e d t e r m s , u s e d a s a m e a s u r e o f t h e s e n s i t i v i t y o f t h i s p r o ce du r e t o t h e i r d e t e r m i n a t i o n , a r e t a k e n f r o m t h e p o l y n o m i a l fi t w i t h t h e li n e a r t e r m ss ub t r a c t e d p r i o r t o r e m o v i n g t h e i n s i g n i fi c a n tt e r m s . T h e v a l u e s g i v e nh e r e f o r t h e n o n - ze r o m u l t i p o l ec o e ffi c i e n t s a r e o b t a i n e d i n t h e fin a l s t e p o f t h e p r o ce du r e , i. e . t h ec o n s t r a i n e dfi t w i t h t h e li n e a r t e r m ss ub t r a c t e d . V e r s i o nS o l v e r t i m e M o d e l s i z e σ R e s i du a l s b b b b ( m i nu t e s × c o r e s )( N r o f n o d e s )( T m )( T m / m )( T m / m )( T m / m )( T m / m ) . × . M . × − ( ) ± . × − . ± . × ( ) ± . × − . ± . × . × . M . × − ( ) ± . × − . ± . × ( ) ± . . ± . × . × . M . × − ( ) ± . × − . ± . × ( ) ± . × . ± . × . × . M . × − ( ) ± . × − − . ± . × ( ) ± . × . ± . × . × . M . × − ( ) ± . × − ( ) ± . × − ( ) ± . × . ± . × . × . M . × − ( ) ± . × − − . ± . ( ) ± . × . ± . × .3.2. Deviations from Maxwell’s Equations Figure 19 shows the evolution of the diagnostic criterion developed in Sec. 4.2
Version 1 Quad Type A Map type nodalVersion 2 Quad Type A Map type nodalVersion 3 Quad Type A Map type nodalVersion 4 Quad Type A Map type nodalVersion 5 Quad Type A Map type nodalVersion 6 Quad Type A Map type nodal
X (m) I n ve rs e o f r oo t o f s u mm e d s q u a r e s o f r e l a ti ve e rr o rs i n C u r l Z ( × - ) Figure 19: Improvement in the consistency of the field calculations with Maxwell’s equationsusing the diagnostic criterion developed in Sec. 4.2 for the six finite-element models describedin Sec. 5.2. as a quantitative assessment of the level of numerical deviations from Maxwell’sequations. The relationship of these results to the particular mesh refinementsis discussed above in Sec. 5.2. The refinement on the XZ and 45-degree planeswithin 2.5 mm of the magnet axis (version 2) resulted in a reduction of the Zcomponent of the curl by a factor of more than 20. The introduction of a “micro-beam” surface within 2 mm of the magnet axis of cell size 0.4 mm (version 5)significantly improved the accuracy in the region occupied by the positron beam,which is of the order of 1 mm wide wide in the south arc of CESR and is focusedto a point near the end of each quadrupole. The final step (version 6) introduceda fine mesh volume in the fringe field region, significantly reducing the deviationsfrom Maxwell’s equations in the fringe region, but showing only a small effecton the field integrals.
6. Conclusions and Outlook
The methods developed during this work provide useful input to two types ofaccelerator optics design, implementation, and operational diagnostic projects:29 quantitative assessment of the accuracy of field maps derived from finite-element models provides a diagnostic tool for the use of these field maps inparticle tracking methods for determining accelerator element transportmatrices. These can prove useful at any stage of accelerator developmentfrom calculating the effects of approximations and errors such as misalign-ments during the design process to the diagnosis of deviations in measuredoptical functions from design. • the reliable derivation of the multipole content in discrete field maps en-ables the implementation of higher multipole coefficients in elements oflattice models for accelerator optics. The orders-of-magnitude improve-ment in computing time makes many types of design and diagnostic stud-ies possible when field-map tracking is impractical. We have presenteda case of multipole determination for a quadrupole magnet design withvery small multipole content. The methods described here are generallyapplicable to finite-element models for designs with large multipole termsarising from, for example, severe space constraints. In such cases, boththe tracking and multipole analysis will be of importance for both designand operational optics corrections.A prime example of the need for such studies are the splitter/combiner beam-lines for the Cornell Brookhaven Energy-recovery-linac Test Accelerator [16, 17].The compact nature of these beam-lines, together strict transverse uniformityspecifications, resulted in the use of 3D tracking models for the 42, 78, 114and 150 MeV electrons for the magnet designs. In fact, the design criteria wereachieved through the introduction of angled chamfers on the ends of the poles forthe purpose of shaping the fringe fields, based on the tracking results [18, 19]. Inaddition to the two element model improvements cited above, Bmad provides forfringe field shape parameters which were used effectively in the CBETA latticemodel .An extensive study comparing various methods of producing and applyingtransfer maps, including from discrete field tables, has been performed for thenon-scaling fixed-field alternating-gradient accelerator EMMA [20]. Much workon avoiding numerical noise in generating reliable transfer maps from discretefield maps by using Maxwell’s equations to interpolate inward from a boundingsurface far from the magnet axis has also been done [21, 22, 23, 24, 25], andcitations therein.We look forward in the near future to applying the model refinement anddiagnostics described here to both field-map tracking and implementation ofmultipole content in the magnetic lattice elements for the relatively recentlyinstalled south arc region of CESR. While we do not have specific expectationsof surprises, we find value in the results of such a first-time study, even whenthey affirm the validity of past approximations for some applications. Private communication from J.S. Berg of Brookhaven National Laboratory . Acknowledgments We wish to acknowledge useful conversations with I. Bazarov, G. Hoffstaet-ter, V. Khachatryan, D. Rubin, D. Sagan and S. Wang, of CLASSE, as well aswith J.S. Berg, F. M´eot and N. Tsoupas of Brookhaven National Laboratory.Valuable assistance was also provided by Y. Zhilichev of Dassault Syst`emes,including a critical reading of Sec. 3. J. Shanks provided enlightenment con-cerning CHESSU optics design. Help was also provided by N. Verboncoeur inthe context of the National Science Foundations Research Experience for Under-graduates program. This work is supported by the National Science Foundation.
References [1] G. Strang and G. Fix,
An Analysis of the Finite Element Method, SecondEdition (Wellesley-Cambridge, 2008).[2] O. Zienkiewicz, R. Taylor, and J. Zhu, eds.,
The Finite Element Method:its Basis and Fundamentals (Seventh Edition) , seventh edition ed. (ElsevierButterworth-Heinemann, Oxford, 2013).[3] J. Akin,
Finite Element Analysis with Error Estimators (Butterworth-Heinemann, Oxford, 2005).[4] J. Shanks, J. Barley, S. Barrett, M. Billing, G. Codner, Y. Li, X. Liu,A. Lyndaker, D. Rice, N. Rider, D. L. Rubin, A. Temnykh, and S. T.Wang, Phys. Rev. Accel. Beams , 021602 (2019).[5] K. Berkelman, A Personal History of CESR and CLEO, The Cornell Elec-tron Storage Ring and Its Main Particle Detector Facility (World Scientific,New York, 2004).[6] D. Sagan, in
Proceedings of the 8th International Conference on Compu-tational Accelerator Physics, ICAP 2004, St. Petersburg, Russia, June 29-July 2, 2004 , Vol. A558 (2006) pp. 356–359.[7] F. M´eot, J. Berg, S. Brooks, D. Trbojevic, N. Tsoupas, and J. Crittenden,Int. J. Mod. Phys. A , 1942038 (2019).[8] F. M´eot, S. Brooks, J. Crittenden, D. Trbojevic, and N. Tsoupas, in Proc. 13th International Computational Accelerator Physics Conference(ICAP’18), Key West, FL, USA, 20-24 October 2018 , International Com-putational Accelerator Physics Conference No. 13 (JACoW Publishing,Geneva, Switzerland, 2019) pp. 186–190.[9] J. Crittenden, V. Khachatryan, and S. Wang, “Magnetic field calculationsfor the CHESS-U combined-function magnet,” CLASSE (2020).[10] J. Simkin and C. W. Trowbridge, in
Proceedings of the 1976 CompumagConference, Oxford, United Kingdom (The International Compumag Soci-ety, 1976) pp. 5–14. 3111] F. James and M. Roos, Computer Physics Communications , 343 (1975).[12] F. James, “ MINUIT , function minimization and error analysis,” CERNprogram library writeup D506 (1998).[13] F. James, “The interpretation of errors,” CERN Geneva, Switzerland(2004).[14] R. Brun, O. Couet, C. E. Vandoni, and P. Zanarini, Computer PhysicsCommunications , 432 (1989).[15] R. Brun, O. Couet, C. Vandoni, P. Zananni, and M. Goossens, “PAW++,Physics Analysis Workstation User’s Guide,” CERN program library longwriteup Q121 (1999).[16] A. Bartnik et al. , Phys. Rev. Lett. , 044803 (2020).[17] C. Gulliford et al. , (2020), accepted for publication in Phys. Rev. ST Accel.Beams, arXiv:2010.14983 [physics.acc-ph] .[18] J. Crittenden et al. , in IPAC2018: Proceedings of the 9th International Par-ticle Accelerator Conference, Vancouver, BC, Canada , Paper THPAF019(2018).[19] J. Crittenden, D. Burke, Y. Fuentes, C. Mayes, and K. Smolenski, in (2017) p. MOPOB59.[20] Y. Giboudot and A. Wolski, Phys. Rev. ST Accel. Beams , 044001 (2012).[21] C. Mitchell and A. J. Dragt, in Proc. of International Computational Ac-celerator Physics Conference (ICAP’15), Shanghai, China, 12-16 October2015 , International Computational Accelerator Physics Conference No. 12(JACoW Publishing, Geneva, Switzerland, 2016) pp. 142–146.[22] C. Mitchell and A. Dragt, in