Helmholtz线圈磁场的仿真研究

1、一、仿真题目要求为了获得一定区域上的匀强磁场，可采用多组Helmholtz线圈结构。一种两对线圈的结构如下列图所示。线圈半径a1，a2，线圈间距h1，h2，以及线圈中通过电流i1，i2可变化量，如图1a所示。为了定量衡量关注区域的磁场均压程度，过轴线做截面ABo1o2，取CD=0.8×AB和Eo3 = 0.8×o1o2，在CD和Eo3线段上每边均匀取20采样点，从而形成如图1b所示的采样节点，定义z方向B的不均压系数为：其中，为所有采样点的z方向磁感应强度平均值；为第n个采样点的z方向磁感应强度值，N为采样点总数。如给定a1 = 1 m， i1=1 A，且h1 > 0

2、.2 m， 情况下，如果要使得到达最小，线圈半径a2，线圈间距h1，h2，以及线圈中通过电流i2应如何取值？(a)线圈结构示意图b磁场采样节点示意图二、公式推导1.对于单线圈情况推导磁感应强度由于仿真要求在xoz平面上取点，对磁感应强度z轴分量进行不均压度分析，那么只需关注在xoz平面上点的，推导过程如下：如上图所示，位于处的电流元，对场点P(r,z)处产生的磁感应强度z轴分量那么整个线圈电流元中从0到，进行积分得到整个线圈对场点P产生的磁感应强度z轴分量2.不均匀度公式的简化题目中给定的不均匀度为所有采样点的磁感应强度标准差除以平均值的绝对值，而由数学推导可得如下化简形式，即采样点磁感应强度

3、平方的均值除以均值的平方再减一，即得不均匀度的平方。此处使用不均匀度的平方是因为可以消去绝对值运算和开方运算，减少运算次数、同时也使目标函数变得可导。三、线圈分布规律猜想根据之前的推导易于得出，对一个中心位于柱坐标系中z=h处、半径为R、通电流为I的线圈，对于轴线上、高度为z处的场点，磁感应强度z向分量为由于该函数是关于z的偶函数，奇数阶导数一定都恒等于0，对其求偶数阶导数如下。而对于4个线圈，分别按照题目中定义的参数如下表:线圈1234中心位置-h2-h1h1h2半径R2R1R1R2电流I2I1I1I2可以求得4线圈同时作用时，原点处磁感应强度的z向分量各阶导数为 而对同一个上述线圈，仍根据

4、之前得出的式子，可以推导出对于与轴线相距x、高度为0处的场点，磁感应强度z向分量为由于问题本身的物理背景，仍然可以看出此函数是关于x的偶函数，奇数阶导数一定都恒等于0。并且根据高等数学的知识，可以证明此函数对于x的各阶导数相当于积分内的函数先对x求偏导，再做对的积分。求导工作量很大，借助MATLAB的符号工具箱对其求原点(0,0)处偶数阶导数的值，如下所示。仍与刚刚一样代入4个线圈的参数，得到4个线圈共同作用时原点处的各阶导数如下 至此，得到了(1)(4)式。令它们均等于0，可以联立求解出4个未知变量R2、h1、h2、I2。联立的方程组如下。得到一组解：R2 = 2.1535923785202

5、54h1 = 0.293332954474377h2 = 1.967170747875403I2 = 3.000000000000000可以认为按照此规律分布的磁场比拟均匀，故可以用此作为优化初值。四、仿真程序编写1.圆周上小段直电流产生的磁场函数dBz(,R,r)在柱坐标系中，按照上述推导，编写了在(R,)处长为Rd的一段电流，方向为，对于位于(r,z)处的场点，产生的磁场大小在z轴上分量的计算函数如下没有乘系数/4。dBz = (theta) (R-r.*cos(theta)./sqrt(r.2-2.*r.*R.*cos(theta)+R.2+z.2).3)2.一个圆产生的磁感应强度函数B

6、z_single(r,z,R,I)在柱坐标系中对位于(r,z)处的场点，在原点处的半径为R的线圈产生的磁场z轴分量计算函数如下。由于每次计算磁场都会乘系数/4，这局部可以留到函数外面计算，故函数的返回值实际上是没有乘系数/4的磁场强度。函数采用自己编写的梯形公式进行积分，其中step为积分步长，last为上一积分步的被积函数值，now为当前积分步的被积函数值。function Bz = Bz_single(r,z,R,I)% to calculate Bz*4*pi/mu of point (x,z) in magnetic field,% with integral steps vector

7、 thetadBz = (theta) (R-r.*cos(theta)./sqrt(r.2-2.*r.*R.*cos(theta)+R.2+z.2).3);step = pi/20;Bz = 0;last=dBz(0); % to store dBz of last stepfor theta = step:step:2*pi; now=dBz(theta); % to store dBz of this step Bz=Bz+(now+last); last=now;Bz=I*R*0.5*step*Bz;end对此函数进行正确性测试。在?大学物理学(A版)·电磁学?中有一道习题1

8、7.6，求一个半圆形电流I在半圆轴线上距离圆心z处的磁场强度。习题答案中给出的磁场强度轴线方向分量为，为了适应此函数，也即。由于题目中是半圆，相当于积分只从0积分到，在不同的step下得到的积分结果如下表step/1000/100/10准确值积分结果z=03.141592 653589 7933.141592 653589 7943.141592 653589 7933.141592 653589 793z=11.110720 734539 6001.110720 734539 5911.110720 734539 5921.110720 734539 592z=20.280992 58924

9、1 6300.280992 589241 6300.280992 589241 6290.280992 589241 629z=30.099345 882657 9590.099345 882657 9610.099345 882657 9610.099345 882657 961可以看出首先函数编写的是准确的，但在这种情况下，积分步长对函数值变化的影响并不是很大，甚至出现了积分步长越粗大，结果反而更精确的反常现象。这有可能是因为除以一个大数之后变得不精确的原因。为了找到一个经济适宜的积分步长，之后再取一个更有一般性的例子，如线圈半径R=1，电流I=2，场点距轴线距离r=0.3，场点距圆心高度

10、z=0.5，这个例子中的各个参数和最终要解决的问题是一个数量级的。测试不同步长情况下的原函数从0积分到2计算一个线圈磁感应强度的积分结果如下表step/10000/1000/100/20/10/5积分结果8.957661 590821 1028.957661 590821 1028.957661 590821 0998.957661 590821 0998.957661 590831 1958.957675 492861 105误差<1e-15<1e-151e-141e-141e-111e-4此处对的步长选取相当于是将圆等效成了多边形，可以看出/20、/10的步长已经很能满足一般的

11、计算精度要求了，也很显然可以从直观上感觉，40边形和20边形对于圆的近似程度是很高的。因此最后的程序中采用步长/20进行计算。再对函数进行效率优化。由于此积分函数在最终的程序中会被反复调用，此函数耗时一定要短。最初想用MATLAB自带的高斯-拉布拉托积分quad()或quadl()，但由于其中涉及高阶差分和自适应步长等过程，精度虽高，却不如简单的梯形积分快。且由之前过程可看出，只要步长适宜，梯形积分也可以满足精度要求。此函数中在自己编写梯形积分时，每个对应的被积函数仅计算一遍，相同的乘法局部均提出到for循环外，没有任何冗余运算，同时也在保证效率前提下将占用变量即系统内存尽可能的进行了缩减。3

12、.两对圆电流总磁感应强度函数Bz_total(R2,h1,h2,I2,r,z)其中R2为原图中外侧线圈半径，h1、h2分别为内、外侧线圈到原点的距离，I2为外侧线圈电流，以和I1同向为正，r、z分别为采样点的半径坐标向量、纵坐标向量。本函数4次调用单个圆电流产生磁场的函数，在柱坐标系中将4个线圈的总效应相加，调用一次即可计算所有采样点的z向磁感应强度，结果返回给矩阵变量Bz。编写的程序如下。function Bz = Bz_total(R2,h1,h2,I2,r,z)R1 = 1; I1 = 1;% if use parallel computingX,Z=meshgrid(r,z);Z1 =

13、 Z+h2; % coil 1Z2 = Z+h1; % coil 2Z3 = Z-h1; % coil 3Z4 = Z-h2; % coil 4Bz1 = Bz_single(X,Z1,R2,I2);Bz2 = Bz_single(X,Z2,R1,I1);Bz3 = Bz_single(X,Z3,R1,I1);Bz4 = Bz_single(X,Z4,R2,I2); Bz = Bz1+Bz2+Bz3+Bz4;end 值得一提的是，本函数没有采用大家普遍采用的for循环方式来计算400个采样点，而是采用了meshgrid函数生成400个点的坐标矩阵网格，然后对坐标矩阵进行操作，并行的同时计算40

14、0个采样点。在不利用其他并行工具的前提下，MATLAB中for循环是一种低效率的串行方式，而通过上网查找资料发现较高版本的MATLAB中矩阵三角函数运算、矩阵开方运算、矩阵点乘除乘方运算都是自动将运算量分配到多个计算单元上并行执行的，能够显著提速。由于我们用到的运算也就是以上几种的组合，因此本函数中先生成了关于r、z坐标向量的两个二维网格矩阵，再直接对这两个网格矩阵进行运算，以到达并行计算的目的。经过实测，其他条件相同的情况下，采用for循环方式做一次优化至3分钟后仍未出结果其他组同学的程序平均时间一般也在10分钟左右；采用坐标网格矩阵的并行方式做一次优化只需要2秒多。4.不均匀度函数diff

15、er(x)由于题中所给采样区的磁场分布一定是对称的，故此处不均匀度函数仅选取了题目中一半仿真区域11×21个点，半径方向20等分取21点，纵轴方向10等分取11点。函数返回值d为所有采样点不均匀度的平方这对于找其最小值并无影响，输入变量x为一个4维列向量，分别是4个待定参数R2、h1、h2、I2。按照之前推导的简化计算量的公式，编写程序如下。function d = differ(x)% x=R2;h1;h2;I2% d 为不均匀系数delta的平方R2=x(1);h1=x(2);h2=x(3);I2=x(4);% 每边采样点数设置r = 0:0.04:0.8;z = 0:0.08*

16、h1:0.8*h1; % z向仅对上半局部进行采样Bz_half = Bz_total(R2,h1,h2,I2,r,z);Bz = Bz_half;Bz_half(2:end,:);d = mean2(Bz.2)/mean2(Bz)2-1;end因为此函数中没有绝对值，之前的其他函数局部也没有超越函数，所以至此，待优化的目标函数是可导的，有利于后续优化算法的选择。5.优化主函数选用优化工具箱中的函数fmincon进行优化，找到目标函数的最小值。其中约束条件的设置在矩阵A1中和向量v1、v2中，限定了h2>h1、h2>0.2、-3<I2<3。向量v1、v2的其他局部为人为

17、设定的另外约束，为了方便找到最优解。Clear all; clc;tic;mu=4*pi*1e-7;R1=1;I1=1;% variable x=R2;h1;h2;I2;x0=1;0.5;1;3;A1=0,1,-1,0;b1=0;v1=0.1;0.2;0.3;-3;v2=4;0.7;4;3;opt=optimset(largescale,off,MaxIter,200,MaxFun,1000);x,fv,ef,out,grad,hess=fmincon(differ,x0,A1,b1,v1,v2,opt)R2=x(1),h1=x(2),h2=x(3),I2=x(4)delta=sqrt(fv)

18、toc; 五、数值试验1.按之前导数为零的推导选取优化初值首先选取了之前按照原点处各阶导数为0推导出的一组参数作为初值，变量范围设在其附近，得到了一组优化结果，画出其磁感应强度分布如下列图(a)。但随后发现当初值取为其他值时会出现不均匀度更小的解，如下列图(b)。这说明之前的推导是不正确的，它仅仅保证了原点附近磁场足够均匀，却不能保证整个采样区域的不均匀度最小。(a)按照导数为零得到的结果(b)按照另一组初值得到的结果2.按数值试验逐渐缩小范围的方法选取初值我们希望能通过数值实验先计算几组参数值下的不均匀系数平方，再进一步缩小范围寻找其最小值，确定所优化目标的大致位置。由于要优化四个参数，假设参数值一起改变，不太方便观察，可先固定两个参数不变，改变其他两个，再反过来改变另两个即可。根据之前几次优化发现，一般最后结果都会有h1=0.2，I2=3，而另外两个量