Example#6: Simple Quadrupole Magnet

Version : 0.2

Introduction

This example considers a simple quadrupole. Note that it has been derived very quickly and its design has not been optimized at all. This example is in some way similar to the simple dipole magnet of Example#5.nb and all the remaks made in that case are valid here. It is an iron dominated structure. The user should at least first read the introduction of Example#5.nb before playing with this one. The quadrupole presents hyperbolic pole faces with a chamfer.

 

Load and Initialize Radia

The following instruction loads the Radia package and returns the Radia version number. It also turns off some warnings.

 

<<Radia`;
Off[General::"spell1"];

 

Defining The Model

A routine to create the geometry.

 

Nn=4;
ironcolor={0,0.5,1};coilcolor={1,0,0};
depth=1.2*width/2;
If[hyp<0.2,hyp=0.2];

geo[]:=(
rap=0.5;
ct={0,0,0};
z0=gap/2;y0=width/2;
amax=N[hyp*ArcSinh[y0/z0]];dz=z0*(Cosh[amax]-1);
qq=N[Table[{z0*Sinh[a/hyp],z0*Cosh[a]},{a,0,amax*(1+2/np),amax/np}]];
hh=qq[[np+1]][[2]]+height*rap-dz;
qq[[np+2]]={qq[[np+1]][[1]],hh};
qq[[np+3]]={0,hh};
g1=radObjThckPgn[thick/4,thick/2,qq];
 radObjDivMag[g1,n1];
 
g2=radObjRecMag[{thick/4,width/4,gap/2+height*(1/2+rap/2)},{thick/2,width/2,height*(1-rap)}];
radObjDivMag[g2,n2];

gg=radObjCnt[{g1,g2}];  
gp=radObjCutMag[gg,{thick/2-chamfer-gap/2,0,0},{1,0,-1}][[1]];

g3=radObjRecMag[{thick/4,width/4,gap/2+height+depth/2},{thick/2,width/2,
          depth}];

cy={"cyl",{{0,width/2,gap/2+height},{1,0,0}},{0,0,gap/2+height},
        2*depth/width};
radObjDivMag[g3,{nr3,np3,nx},cy];

tan=N[Tan[2*Pi/2/Nn]];
len=tan*(height+gap/2) -width/2;
g4=radObjRecMag[{thick/4,width/2+len/2,gap/2+height+depth/2},{thick/2,len,
          depth}];
radObjDivMag[g4,n4];

posy=width/2+len;posz=posy/tan;
g5=radObjThckPgn[thick/4,
        thick/2,{{posy,posz},{posy,posz+depth},{posy+depth*tan,posz+depth}
 }];
cyl={"cyl",{{0,posy,posz},{1,0,0}},{0,posy,posz+depth},1};
   
radObjDivMag[g5,{nr5,np5,nx},cyl];

Rmax=Rmin-width/2+gap/2+offset-2;
coil1=radObjRaceTrk[{0,0,gap/2+height/2+offset/2},{Rmin,Rmax},{thick,width-2*Rmin},
        height-offset,3,CurDens];
 radObjDrwAtr[coil1,coilcolor];

hh=(height-offset)/2;

coil2=radObjRaceTrk[{0,0,gap/2+height-hh/2},{Rmax,Rmax+hh*0.8},{thick,width-2*Rmin},
        hh,3,CurDens];
 radObjDrwAtr[coil2,coilcolor];
 
g=radObjCnt[{gp,g3,g4,g5}];
radObjDrwAtr[g,ironcolor];
  
RadTrfZerPerp[g,ct,{1,0,0}];
RadTrfZerPerp[g,ct,{0,1,0}];
t=radObjCnt[{g,coil1,coil2}];
RadTrfZerPara[t,ct,{0,N[Cos[Pi/Nn]],N[Sin[Pi/Nn]]}];

radTrfMlt[t,radTrfR
ot[ct,{1,0,0},4*Pi/Nn],Nn/2];
radMatApl[g,ironmat];
 radTrfOrnt[t,radTrfRot[{0,0,0},{1,0,0},N[Pi/Nn]]];)
 
(* ------  For the Harmonic Analysis --------- *)
 
harm[Obj_,{y_,z_,ro_},Np_]:=Module[{},{Np,ro,Table[
radFldInt[Obj,"inf","ibz",N[{-1,y+ro*Cos[tet],z+ro*Sin[tet]}],N[{1,y+ro*Cos[tet],z+ro*Sin[tet]}]]+
radFldInt[Obj,"inf","iby",{-1,y+ro*Cos[tet],z+ro*Sin[tet]},{1,y+ro*Cos[tet],z+ro*Sin[tet]}]*I,{tet,0,2*Pi*(1-1/Np),2*Pi/Np}]}]

multipole[w_List,q_]:=( w[[3]].Table[Exp[-2*I*Pi/w[[1]]*q*p],{p,0,w[[1]]-1}]/w[[1]]/w[[2]]^q)

 

Creating the Model and Solving

Creating and solving the magnetic structure. Note that the results are strongly dependent on the segmentation parameters. nx and ny are easy to use numbers which determine the segmentation in the length of the qaudrupole and in the plane of the pole. Change them and watch the periodic part of the iron to see what segmentation it represents. Note it may happen (in rare cases) for some geometry and values of nx and ny that the relaxation is either unstable or uncomplete. This will be reported in the result of the computation, the cure is typically to reduce one of the index nx or ny and increase the other one.
The shape of the pole is a hyperbole approximated by 2 x np flat segments. The parameter np is important if one cares of removing the dodecapole over a large range around the axis. The parameter hyp define the shape of the pole faces. hyp=1 is a perfect hyperbolic shape, lower values result in flat poles, higher values gives poles with more curvature. One easily check that even with a low segmentation the dodecapole is minimzed for a value of hyp close to 1 as expected for a non saturated quadrupole.
Taking as a reference a computation made with nx=10; ny=6; n1={nx,ny,{10,1}}; n2=n1; np3=4; nr3=ny; n4={nx,8,ny}; np5=4; nr5=ny; occupying 190 Mbytes and taking 6 hours of cpu time one obtains a gradient of 15.43 T / m. The segmentation parameters used below gives the gradient with a 1% precision within 17 seconds of solving and little memory.

 

np=4;                  (* see above *)
hyp=1;     (* must be > 0.2 *)

gap=40;     (* magnetic gap *)
width=30;    (* pole width *)
height=50;    (* height of iron in coil *)
thick=60;    (* length of Quad *)
chamfer=8;    

(* coils *)
CurDens=-3;    (* current density A/mm2 *)
Rmin=2;offset=10;

(* segmentation parameters *)

nx=2;ny=2;
n1={nx,ny,{3,2}};
n2={nx,ny,{3,1}};
np3=2;nr3=ny;
n4={nx,3,ny};
np5=2;nr5=ny;

t0=AbsoluteTime[];
radUtiDelAll[];
ironmat=radMatSatIso[{2000,2},{0.1,2},{0.1,2}];
geo[];
size=radObjDegFre[t];

t1=AbsoluteTime[];
Nmax=1000;
res=RadSolve[t,0.00001,Nmax];
t2=AbsoluteTime[];
Bz=radFld[t,"Bz",{0,1,0}]*1000; 
Iz=radFldInt[t,"inf","ibz",{-1,1,0},{1,1,0}];
Iz1=radFldInt[t,"inf","ibz",{-1,10,0},{1,10,0}]/10;

Print["Mag_Max  H_Max  N_Iter = ",N[res[[2]],3]," T ",N[res[[3]],3]," T ",res[[4]]];
If[res[[4]]==Nmax,Print["Unstable or Uncomplete Relaxation"],];
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[" Gradient = ", N[Bz,4]," T/m "];
Print["Int. Quad. @ 1 mm = ", N[Iz,5],"T"];
Print["Delta Int. Quad. @ 10 mm = ", N[100*(Iz1/Iz-1),2],"%"];

 

Harmonic Analysis

Sequence of Integrated Dipole [T mm], Quadrupole[T], Sextupole [T / mm],.. The real (imaginary) part is the normal (skew) component. The functions harm and multipole are dedined in the section "Defining The Model". As expected, the next dominant multipole following the quadrupole is the dodecapole.

 

nharm=10;radius=2;
w=harm[t,{0,0,radius},nharm];
N[Table[multipole[w,i],{i,0,nharm-1}],4] 

 

QuickDraw 3D Graphics

Here are instructions to plot using QuickDraw3D the whole quadrupole or the periodic part of the iron structure.

 

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

(* The Periodic part of the Iron  *)
radObjDrwQD3D[g]; 

 

Mathematica Graphics

The plot of the whole quadrupole requires significant memory for the Kernel. Depending on your memory you may or may not be capable of plotting it. However, one may just want to plot the periodic part of the iron structure to see its segmentation. To plot the whole quadrupole replace "radObjDrw[g]" by "radObjDrw[t]" and be preapred to wait for some time. This time can be minimized if you only plot loosely segmented structures.

 

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

Magnetic Field Plots

Here are a few plots of magnetic field and field integrals.

 

ob=t;z=0;x=0;ymax=40;Npt=20;
RadPlotOptions[];
bz=radFldLst[ob,"Bz",{x,0,z},{x,ymax,z},Npt,"arg",0];
ListPlot[bz,PlotJoined->True,
FrameLabel->{"Y [mm]","Bz [T]","X = 0 mm, Z = 0 mm",""}]; 

 

ob=t;z=0;ymax=10;
RadPlotOptions[];
Iz0=radFldInt[ob,"inf","ibz",{-1,1,z},{1,1,z}];
Plot[(radFldInt[ob,"inf","ibz",{-1,y,z},{1,y,z}]/(y*Iz0)-1)*100,{y,0,ymax}
,FrameLabel->{"Y [mm]","dIz [%]"," Z = 0 mm",""}];