Example#5: Simple Dipole Magnet

Version : 0.3

Introduction

This example concerns the geometries dominated by iron. At any point in space one can consider that the magnetic field is the sum of two contributions. One comes from the real sources of magnetic field such as coils or permanent magnets. The other comes from the iron through its magnetization. Even though the magnetization of the iron is induced by the coils or permanent magnets, its field contribution to the point of interest can be much larger than the direct contribution form coils and magnets. We call these geometries: Iron Dominated. Electro-magnet dipole and quadrupole used in particle accelerators fall into this category. Most of undulators and wigglers used as synchrotron radiation sources do not fall into this category. The field computation with Radia in the case of Iron Dominated geometries presents a few specific difficulties and is usually less accurate than in the case of coil or permanent magnet dominated structure. Nevertheless, special methods have been developed to reach a reasonable precision within a reasonable cpu time and memory. The example below is that of a simple dipole steerer made of a closed circuit of iron with a small gap and a coil wounded around the circuit that drives some flux in the iron (see the graphics below). This example is more delicate than all previous examples and we advise the beginner to first get experienced with them before diving into this one.

Among the most important things to remember if one wants to reach a good precision in a reasonable time are:

- Always segment the corners of the iron circuits as parallel or perpendicular to the flux line as possible. For right angle corners, this can be done with the circular or ellipsoidal mode of segmentation (see below). This is extremely important. In this example, we make use of the circular segmentation twice (the other two corners are segmented similarly by symmetry) it is made in the section entitled "Defining the Model" when calling "radObjDivMag[g3,{nbr,nbp,n3[[1]]},typ]" and "radObjDivMag[g5,{nbr,nbp,n5[[1]]},typ]".

- Use a narrower segmentation on the iron faces close to the region of interest.

- Start with low segmentation numbers and increase them gradually to check that the field is stable. Beware that the memory and cpu time tends to grow like the square number of elements.

- Use symmetries as much as possible. It saves both memory and CPU time. The steerer dipole shown below has a symmetry of order 2 x 2. 16 times more memory and 4 times more cpu time would have been needed without using the symmetries.

As for the previous examples, try to modify some parameters (see the Reference Guide for the explanations on all Radia functions) and re-execute the corresponding section. All sections must be evaluated in the order of presentation. A section may be evaluated several times with the same or different parameters if the previous sections have been evaluated before. The only limitation is memory. The Radia.exe memory can be re-initialized by re-executing the section entitled "Load and Initialize Radia". The Kernel memory can be re-initialized by executing Exit[] and then re-executing all sections of this example. If the Front-End is running out of memory, close some windows or close the Front-End and start it again. Note that if by mistake any section is executed before the previous ones are executed, it may be necessary to exit the Kernel and re-start everything from the very beginning. The maximum precision requires large memory and long cpu time and ultimately, you may want to add extra memory to your computer. However remember that the memory required tends to be proportional to the square number of sub-elements and increasing the segmentation by 2 in all three directions results in 8 times more sub-elements and 64 times more memory. Consequently, one should increase the segmentation number in some "intelligent" manner, this requires a bit of experience. One way to develop this skill is to watch the precision following an increase of the segmentation in one object and one plane at a time. By this method, one quickly identifies the objects and the directions which require more segmentation to improve the accuracy.

Load and Initialize Radia

The following instruction loads the Radia package and returns the Radia version number.

 

<<Radia`;

Defining The Model

A routine to create the geometry.

 

eps=0;     (* used for debugging *)
Off[General::"spell1"];
geom[circ_]:=(
ironcolor={0,0.5,1};
coilcolor={1,0,0};

(* The pole faces *)

lx1=thick/2;ly1=width;lz1=20;l1={lx1,ly1,lz1};
k1={{thick/4-chamfer/2,0,gap/2},{thick/2-chamfer,ly1-2*chamfer}};
k2={{thick/4,0,gap/2+chamfer},{thick/2,ly1}};
k3={{thick/4,0,gap/2+lz1},{thick/2,ly1}};
g1=radObjMltExtRtg[{k1,k2,k3}];
radObjDivMag[g1,n1];

(* Vertical segment on top of pole faces *)

lx2=thick/2;ly2=ly1;lz2=30;l2={lx2,ly2,lz2};
p2={thick/4,0,lz1+gap/2+lz2/2+1*eps};
g2=radObjRecMag[p2,l2];radObjDivMag[g2,n2];

(* Corner *)

lx3=thick/2;ly3=ly2;lz3=ly2*1.25;l3={lx3,ly3,lz3};
p3={thick/4,0,lz1+gap/2+lz2+lz3/2+2*eps};
g3=radObjRecMag[p3,l3];
typ={"cyl",{p3+{0,+ly3/2,-lz3/2},{1,0,0}},p3+{0,-ly3/2,-lz3/2},lz3/ly3};

If[circ==1,
radObjDivMag[g3,{nbr,nbp,n3[[1]]},typ],
radObjDivMag[g3,n3]];

(* Horizontal segment between the corners *)

lx4=thick/2;ly4=80;lz4=lz3;l4={lx4,ly4,lz4};p4={thick/4,ly3/2+eps+ly4/2,p3[[3]]};
g4=radObjRecMag[p4,l4];radObjDivMag[g4,n4];

(* The other corner  *)

lx5=thick/2;ly5=lz4*1.25;lz5=lz4;l5={lx5,ly5,lz5};
p5={thick/4,p4[[2]]+eps+(ly4+ly5)/2,p4[[3]]};
g5=radObjRecMag[p5,l5];


typ={"cyl",{p5+{0,-ly5/2,-lz5/2},{1,0,0}},p5+{0,ly5/2,-lz5/2},lz5/ly5};

If[circ==1,
radObjDivMag[g5,{nbr,nbp,n5[[1]]},typ],
radObjDivMag[g5,n5]];

(* The vertical segment inside the coil *)

lx6=thick/2;ly6=ly5;lz6=gap/2+lz1+lz2;l6={lx6,ly6,lz6};
p6={thick/4,p5[[2]],p5[[3]]-(lz6+lz5)/2-eps};
g6=radObjRecMag[p6,l6];radObjDivMag[g6,n6];

(* Generation of the coil *)

Rmin=5;Rmax=40;Nseg=4;H=2*lz6-5;
CurDens=current/H/(Rmax-Rmin);
pc={0,p6[[2]],0};
coil=radObjRaceTrk[pc,{Rmin,Rmax},{thick,ly6},H,3,CurDens];
radObjDrwAtr[coil,coilcolor];

(* make container and set the colors *)

g=radObjCnt[{g1,g2,g3,g4,g5,g6}];

radObjDrwAtr[g,ironcolor];
radMatApl[g,ironmat];
t=radObjCnt[{g,coil}];

(* define the symmetries *)

RadTrfZerPerp[g,{0,0,0},{1,0,0}];
RadTrfZerPara[g,{0,0,0},{0,0,1}];
);

 

Creating the Model and Solving

Creating and solving of the field geometry. The Bz field printed is the vertical field in the middle of the gap. In this computation, the subdivision in the iron is low in order to achieve a short CPU time at the cost of a reduced precision. The excitation current is low and the iron has a very large permeability. Therefore there is no saturation in the iron and the Ampere law gives a good estimate of the field on axis to which the computed field can be compared. The sources of discrepancy between the computed and estimated fields from Ampere's law are the residual small deviation from uniformity in the gap and the unsufficient segmentation. A closer agreement between them is reached if one uses the following higher segmentation parameters : nx=3; nbp=3; nbr=3; n1={nx,3,2}; n2={nx,3,3};n4={nx,3,3} n6={nx,3,3}, the cpu time is longer but still reasonable. One can observe some saturation in the iron by either increasing the current in the coil (for example current =-10000) or using a lower saturated magnetization when calling "radMatSatIso" or by using an other material. The computed Mag_Max and H_Max are the maximum value of the magnetization vector and H vector in the iron circuit both expressed in Tesla. If saturation starts somewhere in the circuit, Mag_Max approaches the saturated magnetization of the material (2T in this example as defined in "radMatSatIso") and H_Max deviates from 0. One may want to vary the gap, thick, width and chamfer parameters but one may need to increase or reduce some segmentation parameter accordingly to keep a reasonable precision for a reasonable cpu time. Just try and play by changing some parameters and re-executing ! The only risk is if one sets a too large segmentation parameter resulting in insufficient memory or too long cpu time. If this happens, in the worst case you will have to kill the Mathematica Kernel by one way on an other, if Radia.exe is still alive kill it also (by switching to it and then exiting) then re-execute all sec tions of this notebook starting by "Load and Initialize Radia" etc... One may estimate the memory requirement by looking at the size of the interaction matrix. You may then switch to the next sections to see how to display the 3D graphics of your geometry or make field or field integral plots at whatever point you are interested in.

 

(* Geometry Parameters *)

gap=10;       (* in mm  *)
thick=50;     (* in mm  *)
width=40;     (* in mm  *)
chamfer=8;     (* in mm  *)
current=-2000;    (* in Ampere  *)

(* Segmentation Parameters *)

nx=2;      
nbp=2;nbr=2;     (* for the corners *)
n1={nx,3,2};    (* pole faces *)
n2={nx,2,2};    (* Small Vert. Arm *)
n3={nx,2,2};    
n4={nx,2,2};    (* Horiz.  Arm  *)
n5={nx,2,2};    
n6={nx,2,2};    (* Inside the Coil *)

(* Build the geometry *)

t0=AbsoluteTime[];
radUtiDelAll[];
ironmat=radMatSatIso[{20000,2},{0.1,2},{0.1,2}];
geom[1];
size=radObjDegFre[t];

(* Solve the geometry *)

t1=AbsoluteTime[];
res=RadSolve[t,0.0002,1500];
t2=AbsoluteTime[];

(* Print The results *)

b0=radFld[t,"Bz",{0,0,0}];
bampere=(-4*Pi*current/gap)/10000;
r=b0/bampere;

Print["Mag_Max  H_Max  N_Iter = ",N[res[[2]],3]," T ",N[res[[3]],3]," T ",res[[4]]];
Print["Built & Solved in ",Round[t1-t0]," & ",Round[t2-t1], " seconds"]
Print["Interaction Matrix : ",size," X ",size," or ", N[size*size*4/1000000,2]," Mbytes"]
Print["Bz = ", N[b0,4]," T,   Bz Computed / Bz Ampere Law  = ",N[r,4]]

QuickDraw 3D Graphics

If you have QuickDraw 3D installed on your platform, execute one of the following lines to display the geometry in a dedicated window where you will be able to zoom, rotate etc... When you have finished, Quit from the Menu File of Radia.exe (for Macs) or Exit from the Menu File of Windows (for PC) to release the memory required by QuickDraw3D and continue working in Mathematica. Each of the call will draw a different object or part of the whole assembly. The color may differ in a sub-object if the color has been applied to the container rather than to each sub-object as it is the case in this model. Calling "radObjDrwQD3D" is the recommended method for debugging a newly written geometry. It may happen that the segmentation lines are not drawn correctly by "radObjDrwQD3D", in case of doubt, the user can make a Mathematica plot (see the section entitled "Mathematica Graphics") which is more time consuming but more accurate and suitable for postscript printing. The QuickDraw3D plot can be saved into a file of type "3DMF" using the "Save As" in the file menu. A number of commercial codes can read such files and export them to an other 2D or 3D format.

 

(* The whole geometry *)
radObjDrwQD3D[t]; 

(* The Iron only  *)
radObjDrwQD3D[g]; 

(* One corner *)
radObjDrwQD3D[g3]; 

 

Mathematica Graphics

These instructions plot the geometry of the steerer dipole.

 

draw=radObjDrw[t];
RadPlot3DOptions[];
Show[Graphics3D[draw]
,ViewPoint->{5,-1.5,1}
,PlotRange->All
,AmbientLight->GrayLevel[0.1]]; 

Magnetic Field Plots

Here are a few plots of field and field integrals in the gap. Note that the plot is made by calling "radFldLst " followed by a "ListPlot". The same result can be obtained by calling "radFld" followed by "Plot" but it would be significantly slower because of the large number of calls made to Mathlink . The mathlink protocol used by Mathematica is responsible for the data exchange between the Kernel and Radia.exe processes. It is rather slow (dependening on platforms). It is also much easier to control the number of points used in the graphichs when calling "radFldLst" rather than "radFld". Each graph corresponds to the field created by the whole steerer. It is possible to compute separately the contribution from the coil or from the iron. To do so, one simply needs to replace ob=t by ob=coil or ob=g. It is also possible to compute the contribution from each part of the iron circuit g1,g2...g6. To do so, one would use ob=g1 or ob=g2 etc.. However, one should remember that the symmetries has been assigned to the container g and not to each element g1,g2..g6 therefore the field from each element must be multiplied by 2 x 2 = 4. One should remember that even though the field from the whole steerer is the sum of the field from the coil + 4 times the field from g1+g2+...g6, the magnetization of any element gi depends on the magnetization of all the other elements. It is nevertheless quite interesting to see which part of the magnetized iron circuit contributes the most to the field in the gap and what is its field gradient.

 

ob=t;y=1;z=1;xmax=30;Nn=40;
RadPlotOptions[];
bz=radFldLst[ob,"Bz",{-xmax,y,z},{xmax,y,z},Nn,"arg",-xmax];
ListPlot[ bz,PlotJoined\[Rule]True,
FrameLabel->{"X [mm]","Bz [T]","Y = 1 mm, Z = 1 mm",""}]; 

 

ob=t;z=3;
pl=Plot3D[radFld[ob,"Bz",{x,y,z}],{x,-thick/2,thick/2},{y,-width/2,width/2}
,AxesLabel->{"X [mm]","Y [mm]","Bz [T]"}]; 

 

ob=t;
Plot[radFldInt[ob,"inf","ibz",{-1,y,z},{1,y,z}],{y,-20,20},
FrameLabel->{"Y [mm]","Iz [T mm]","Z = 0",""}]; 

Magnetic Force

One computes the force over the element g1 induced by the whole structure. One should be cautious that the cpu time to compute the force can be very long and one must set the segmentation parameter k to the lowest values consistent with sufficient precision . Taking into account the symmetry with respect to the x axis {1,0,0} would remove the Fx component and double the Fy and Fz components.

 

k={1,1,2};
t1=AbsoluteTime[];
fr=radFldEnrFrc[g1,t,"",k];
t2=AbsoluteTime[];
Print["{Fx, Fy, Fz} = ",N[fr,3]," Newton"];
Print ["cpu time : ",N[t2-t1,3]," seconds"]; 

 

Creating the Model and Solving with Rectangular Segmentation in the Corner

The same model is created but the segmentation in the corner is made with parallelepipedic sub-elements. The precision of the field in the gap is much poorer. It can be improved by using larger values of n3 and n5 such as {nx,4,4} or {nx,6,6}, but the CPU time and memory requirements grow quickly.

 

(* Geometry Parameters *)

gap=10;       (* in mm  *)
thick=50;     (* in mm  *)
width=40;     (* in mm  *)
chamfer=8;     (* in mm  *)
current=-2000;    (* in Ampere  *)

(* Segmentation Parameters *)

nx=2;      
nbp=2;nbr=2;     (* for the corners *)
n1={nx,3,2};    (* pole faces *)
n2={nx,2,2};    (* Small Vert. Arm *)
n3={nx,2,2};    
n4={nx,2,2};    (* Horiz.  Arm  *)
n5={nx,2,2};    
n6={nx,2,2};    (* Inside the Coil *)

(* Build the geometry *)

t0=AbsoluteTime[];
radUtiDelAll[];
ironmat=radMatSatIso[{20000,2},{0.1,2},{0.1,2}];
geom[0];

(* Solve the geometry *)

t1=AbsoluteTime[];
res=RadSolve[t,0.0001,1500];
t2=AbsoluteTime[];

(* Print The results *)

b0=radFld[t,"Bz",{0,0,0}];
bampere=(-4*Pi*current/gap)/10000;
r=b0/bampere;

Print["Mag_Max  H_Max  N_Iter = ",N[res[[2]],3]," T ",N[res[[3]],3]," T ",res[[4]]];
Print["Built & Solved in ",Round[t1-t0]," & ",Round[t2-t1], " seconds"]
Print["Bz = ", N[b0,4]," T,   Bz Computed / Bz Ampere Law  = ",N[r,4]]

draw=radObjDrw[t];
RadPlot3DOptions[];
Show[Graphics3D[draw]
,ViewPoint->{5,-1.5,1}
,PlotRange->All
 ,AmbientLight -> GrayLevel[0.1]];
 
 Print[" Segmenting iron corners into parallellepiped should be avoided. Instead, one should use a circular or ellipsoidal segmentation"]