亥姆霍兹线圈仿真

题目：为了获得一定区域上的匀强磁场，可采用多组Helmholtz线圏结构。一种两对线囲的结构如图1所示。线圏半径g亚、线囲间距力”力2,以及线圈中通过电流触&可变化童，如图1 3 所示。为了定童衡量关注区域的磁场均压程度，过轴线做截面凡旳辺，取n).8X/LP和耳勺二 ().8XAo3在QP和线段上每边均匀取20采样点，从而形成如图1 (b)所示的采样节点， 定义2方向B的不均压系数为：线圈结构示意图其中，R为所有采样点的2方向磁感应强度平均值；B$)为第个釆样点的2方向磁感应ZZ强度值。2为采样点总数。定51参数：1)1=1、/%, b? =o问题：如呆规定ii = - i2，问b2. a, a?如何取值可以使得d最小，即关注区域黴 场“最均匀”。(b)黴场采样节点示意图图1 两对线圈产生匀强磁场示意图仿真要求：1) 写出给定起点、终点、场点坐标，編制空间中一栽流直线段在任意观察点的磁感应强度计 笄程序。2) 写出单个囲环线圈空间任意点磁感应賤度的计算程序，并进行验证。3) 复习Matlab中优化工具箱的使用。二、仿真与分析：(一) 、一載流直线段在任意观察点的磁感应强度1. 理论分析：如图，在直角坐标系内，设坐标原点为AB中点，AB在n轴上，起A(0,0,-L/2),终点B(0,0,L/2)o根据书中例3-1的结论可知，对于通过电流I的直导线AB,任意观察点P(x,y,z) 到AB的距离为R,作观察点P到AB的垂线交于点H,则P点处产生的磁感应强度E为又有B = sinzAPH - sinzBPHc? 4tcRsill zAPH = coszAsill zB PH = cos 乙 BUol4irRcoszA + coszBa2、仿真分析:为了提高计算效率，这卫编程用matlab计算时需用离散的方式来计算磁感应强度：先计笄一小段直导线dl在观测点处产生的磁场强度d瓦再用蛊加的方法，求出整段载流直线段在观测点处的磁场。如：要算一 小段栽流导线丽在P处产生的踐场时，根据毕奥-萨伐尔定律,Luldlx r这里，由于AB是一小段载流导线，可做一个近似运算：TLloI(Jz X BP Llnldz X AP心】/2(飞时+飞时)图2由此，可得到整一段载流导线在观测点P处产生的磁场为:s=fdB根据以上分析，得到计算一载流直线段的matlab程序如附表。现脸证这种算法得到的结呆与理论分析得到之间的误差：以图1-1为例，假设HA=7, HB=3, HP=5,假设电流1 = 1AO则：(1)由理论推导得到的公式计算:4兀 x 10 7 x 1B =4F75C52 + 72+ V52 + 32)=2.6564584 x IO-8而由matlab用益加的方法来计算吋，将AB分成每段长度为0.001的小段来计算和叠加， 笄得的结呆为：B二2.6564577。两者相对谋差低达10一7级别，可见这种算法与理论分析得到 的解析解几乎相同，所以笄法合理。还可以做出这种方法下在一平面上电磁场的分布情况的图像，如下图1-3:图1-3由图可看出，B的方向与电流方向符合右手爍旋定则，箭头的长度代表磁感应孫度的丸 小，可以看到，越靠近栽流线处B越大。结果合理。（二）单个圆环线圈空间任意点離感应强度:K理论推导：对于单囲环线圏所产生的磁场情况，由于此次仿真研究的问题是在平行于圆环平面上的 磁场不均匀程度，即如下图2-1,对磁感应孫度z轴分童比进行不均压度分析，则只需关注 在xoy平面上点的B二，推导过程如下：建立直角坐标系如上图，以圆环圆心为坐标原点，圆环在XOY平柯上，则根据对称性 我们可以得到，取现察点P （%0卫）有以下关系:_ Poldi x r4兀严dl = (Rsinoda, Rcosada, 0)r = (a Rcosa,_Rsina,z)T -1J kdixF= Rsinada Rcosada 0 a Rcosa Rsiiia z得到:zcosaBx =其中，r = vR2 + a2 + z2 2xRcosa 2、仿真与分析：(2-1)小段电流元疊加法：根据(一)中得到的结果，可用分小段盎加的方法来求得一段栽流导线在空间产生的电 磁场情况，此处，所谓的小段电流元盎加法，就是采取这种方法，根据以直代曲的方法，以 等边多边形来代替囲，这样通过多边形的各边产生的电磁场的叠加，即可得到圆形载流线圏 在空间产生电磁场的情况。具体实现程序见附表。(2-2)梯形积分法求解：理论分析已经得到了圆形栽流线圈在空间分布的计笄公式，可在matlab中用梯形积分 的方法对该情况下的踐场的分布。具体编程见附表。【小结】以上所述两种方法郴可得到单个圆形栽流线圏在空间的分布情况。以下探讨这两种 方法的精度和运算速度，以确定后面进行多个线圏的踐场求解时求解方法的选择。在线囲轴线上，线圏轴线上的磁场理论计算较为简便，故对轴线上的磁场情况进分析: 假设线圈半径R二2,线圏电流I二1A,用理论解析法和以上两种数值方法计算线圈轴线 上各点磁场的绪呆进行归纳，如下表所示(数值方法1：即小段电流元盎加法，此处的计算 将圆分为正100边形进冇计笄；数值方法2：用梯形积分法进冇求解)：Z12345理论计笄结呆(10巧2.247941.110720.536200.280990.16093数值方法1计笄结呆(iff)2.247941.110720.536200.280990.16093数值方法2计笄结呆(icr7)2.246461.109990.535840.280810.16083可见，用数值方法2,即用梯形函数积分法进行求解时，在小数点后5为都与理论计算 相同，赭度非常商；而用数值方法1计笄时，由于采用100等边形来近似圆，计算鉛果与理 论结果也比较相近，但还是存在一定谋差，为了减小这种谋差，只要将圆分细一些，就能够 将这种误差控制在妥求范囤之内。但是，用一百等边形来近似圆时，计笄时间已经明显比直 接用积分函数求解的方法长，若将圆再分得细一些，时间会进一步增加。在到进行优化求解 时，必然多次调用此程序，那么计笄时间必然会成倍的增加。综合考虑，后血的多线圈计笄和不均压系数的计算直接采用使用梯形积分法来求解。最后，在进行下一步分析与仿真之前，根据梯形积分法，做出单个圆形栽流线圈的电磁 场沿z轴分布的情况，如下图所示：从图中可看出，在线圏轴线上，磁场賤度在线圏两侧对称分布，在线圈处，磁场賤度最 大，离线圈越远的地方踐场强度越小，符合现实情况。(三) 对不均匀系数的探讨和仿真计算：由题目可知，若要减小不均压系数，则应使得各个观察点的磁场感应賤度均匀分布。从 以上图形和绪呆中我们可以看出，为了进行深入的研究我们可以先从简单入手，对轴线上的 磁感应賤度进冇分析，由于轴线上的磁感应賤度具有代表性，较为容易计算。所以，我们不 妨先研究轴线上的磁感应賤度分布。而之前已经得到单个栽流圆形线圈产生的磁场的空间分布情况，对两个线囲的情况，选用相同的方法，直接进行盎加即可。对两个相距为2h“ R=2, 1=1 A的两个圆环线囲的磁感应强度进行分析。在直角坐标系 内，以两圆环圆心连线中点为坐标原点，XY平血与圆环所在平面平行，这样的话血=黑仿 真得到磁场随町的变化情况即可得到磁场随01的变化情况。仿真得町为不同值时轴线上B 的分布如下：通过图像可看出，两个栽流线圏所产生磁场的分布情况相当于单个栽流线圈所产生磁场 情况的盎加，在轴线上，磁感应賤度最大值在每个圆环中心附近，而两圆环中心连线的中点 处（即原点）的磁感应强度要小于前者。因此，当线囲半径和电流大小不变，而单一改变线圏距离2町时，产生磁场在轴线上的 分布随2叽的变化情况可看作两个单毕曲线波峰的移动。半径不变时，圆环间距离越大，黴 感应强度最大值点越偏移原点，原点处会产生一个波谷；圆环间距离减小，磁感应賤度最大 值向原点移动，两个波毕重叠为一个波峰。而当磁感应孫度分布为两个波舉和一个波谷时, 若增大圆环半径，磁感应强度分布将会变为一个波弹，其效呆与减小圆环间距类似，只是毕 值变小，变化率变小，即分布更加分散。另外，通过仿真发现，R/h不变时，磁感应賤度的总趋势犬致相同，如 21, 2h=l和 R二2, 2h二2的情况，以及R二1, 2h二1.5和R=2, 2h=3的情况。但是，在比例不变的情况下， 逅变大，黴感应賤度分布会向两侧延伸，波舉、波谷的賤度也会减弱。为了使线圏产生的磁场均匀，应通过移动波弹的位豈，使螯体的曲线更“均匀” 一些。 这一点结论，可应用于之后的四个线囲的磁场情况。以下再来探讨如何使四个线圏产生的磁场在空问中的分布更加均匀。四个线圈的情况与 两个线囲的情况相类似，可将四个线圈看作两对线圏，其磁场叠加的情况依旧可看作波毕移 动的情况。值得注意的是，这两对线圈中的电流方向是相反的，所以产生的磁场方向是相反 的，所以相当于把一对波第反倒过来，或者直接进行相减。先固定两组线圏对的半径，通过调节线圏的距离hr h2,亦即调节题中的伤与02,来探 讨在伤与02大概在什么情况下线圏产生的磁场在空间中的分布最均匀。固定两组线圏的半径分别为逅二0.9、坷二4,通过改变hP h2的值进行尝试，最后得到最 均匀的情况如下：这时,hj=2.5, h2=5o通过以上由试探法得到的仿真图像可以看出，场孫分布在z轴上一定范国内基本是均匀 的，说明调整参数得到均匀的磁场具有可行性。为了得到更精确的参敎值，我们以下对最优 解进行求解。现用matlab工具箱对进行最优化求解。优化时，设定a,=l, I. = -I, =1,优化变量 01=血/色、p2=h2/a2. 03=血/勿，考虑实际情况，设程变童下限为0,用遗传算法进行优化求解，得到最优结果为：01 = 0.373、02 = 9.532、角=0.1654。并根据伤、伤、03求得h1=0.373, a2=6.046, h2=57.63优化得到的不均匀系数6= 0.073。为了再进一步验证此时的磁场相对而言是比较均匀的，做出xcy平面的磁场的情况，如 下图：可见，此时的腹场相对而言比较均匀。三、仿真结论:根据以上的讨论分析与仿真的结果，对于本次仿真，可得出以下结论：当 b= h 1=0.37, b2= h2/a2=9.532 时,a2/aj = 6.046 时,不均匀系数6最小，为 6= 0.073。四、总结与反思：1、此次仿真根据老师提供的方向，由浅到深，层层深入：先研究一根载流导线的问题；再 在一根载流导线的基桔上研究圆形栽流线圈的问題；在研究圆形栽流线圈的问题的时候，我 们又先研究一个线圏的情况，再研究两个线囲的情况，最后再分析四个线圏的情况可想 而知，当问题再深入复杂一点的时候，比如研究多个线圏的情况的时候，我们也可以继续深 入研究了。这种研究方法能把一个比较复杂的问题简单化，不只是这道题的解决方法，也是 我们研究其他问题、解决生活实际困难的一个重耍手段；2、本次仿真试验基本原理较为简单，甚至很多公式祁可以直接从课本和课件中直接找到。而难点在于使用matlab对问題进行优化。在編写过程中使用过多种不同的笄法，但是会遇 到求解时问过长或是寻找不到全局最优解等问题，最终保留了鉛果较好的遗传算法；另外， 通过这次仿真也再次考验了我们查找文献的能力。对于本次仿真的问题，一些文献中也有相 似的研究，通过对这些文献的阅读，能够给我们更好的启发。并且，在查找文献资料的过程 中有发现一些更好的算法，但是较为复杂，由于时间关系并没有深入学习，希望能在日后的 课余时问加深对用madab解决优化问题的学习。3、耳归到本题的结果，为了得到不均匀度的最小值，我们最终得到的结果是% = 1片 / a=0373 , b?=/ a2=9.532 时，a2 / aj = 6.046 ,从中可看出 h2 与趣相差 了9倍，血与如相差了 6倍；也就是说，在a尸1时，为了得到0.8仓02 仓“ 坷二0.48大 小面积的均匀磁场，就需要两个半径为6,两个半径为1的线圏通电，线圈距离场点甚至达 到了 57.63,几乎没有实用价值。所以在进行实际工程设计时，要综合考虑这个问题，有时 候，不一定需要得到最优的不均匀度，在允许的范圉内，可能需要降低对均匀度的要求，以 寻求更优的尺寸。程序附录:此函数求解电流元礙场，P, SI, S2分别为场点, 电流元超点终点function B=current_element_B(P,Sl,S2)rl=P-Sl;r2=P-S2;r3=Sl-S2;B=le-7*(cross(rlr3)/(norm(rl)A3 + cross(r2,r3)/(norm(r2)A3)/2;end此函数用于求解载流直线歳场B=0,0,0;P=0,5,0,delta=O 01;for i7:delta:3-deltaB=current_element_B(P,0,0,i,0,0,i+delta)+B; endB此函数通过电流元法求解单一线圏發场，r、h、 p分别为线圈半径，扌居xoy平面距离，场点 function B=circle_l_B(r,h,p)N=100;alpha=0:2 *pi/N: 2 *jpi;B=0,0,0;forn=l:N pO(:,n)=r*cos(alpha(n),r*sin(alpha(n),h, endfor i=l:N-lB=current_element_B(p,pO(: ,i)：pO(j+l)+B, endBhciMTentelementBCPpO/aipOCiP+B； end此函数通过梯形积分法求解单一线圈磁场，r、 z、R分别为场点柱坐标下r、z轴坐标，线圈半 径function Bz=circle 1 _Bz(r,z,R)dBz = (theta)R. *(R-r. *cos(theta)./sqrt(r. A2-2. *t. *R. *cos(theta)+RA2+zA2)A3);step = pi/20,Bz = 0;last=dBz(0);for theta =step: step:21now=dBz(theta);Bz0=(now+last)*0.5*step;lasEow;Bz=Bz+BzO;endBz=le-7*Bz,此函数通过梯形积分法求解四个线圈牍场和delta,其中 betaO=betal, beta?, beta3 function delta = circle4_Bz(betaO) al=l;a2=beta0 (3)*al;hl=betaO(l)*al; h2=beta0(2)*a2;N=20;r=linspace(0,0.8*al,N), z=linspace(-0.8*li 1,0.8*h 1 ,N), R,Z=nieshgrid(r,z);Zl=Z+h2;Z2 = Z+hl;Z3 = Zhl,Z4=Z-h2;Bzl = circle l_Bz(RZ 1 ,a2);Bz2 = circle l_Bz(RZ2,al);Bz3 = circle l_Bz(RZ3,al);Bz4 = circle l_Bz(R,Z4,a2);Bz =Bzl+Bz2+Bz3Bz4, mesh(R,Z,Bz); delta=sqrt(mean2(Bz.A2)/mean2(Bz)A2-l); end