一、仿真题目要求为了获得一定区域上的匀强磁场,可采用多组Helmholtz线圈结构.一种两对线圈的结构如下图所示.线圈半径a1,a2,线圈间距h1,h2,以与线圈中通过电流i1,i2可变化量,如图1〔a〕所示.为了定量衡量关注区域的磁场均压程度,过轴线做截面ABo1o2,取CD=0.8×AB和Eo3 = 0.8×o1o2,在CD和Eo3线段上每边均匀取20采样点,从而形成如图1〔b〕所示的采样节点,定义z方向B的不均压系数为：其中,为所有采样点的z方向磁感应强度平均值；为第n个采样点的z方向磁感应强度值,N为采样点总数.如给定a1 = 1 m, i1=1 A,且h1 > 0.2 m,情况下,如果要使得δ达到最小,线圈半径a2,线圈间距h1,h2,以与线圈中通过电流i2应如何取值？线圈结构示意图〔b〕磁场采样节点示意图二、公式推导1.对于单线圈情况推导磁感应强度由于仿真要求在xoz平面上取点,对磁感应强度z轴分量进行不均压度分析,则只需关注在xoz平面上点的,推导过程如下：如上图所示,位于处的电流元,对场点P处产生的磁感应强度z轴分量则整个线圈电流元中从0到,进行积分得到整个线圈对场点P产生的磁感应强度z轴分量2.不均匀度公式的简化题目中给定的不均匀度为所有采样点的磁感应强度标准差除以平均值的绝对值,而由数学推导可得如下化简形式,即采样点磁感应强度平方的均值除以均值的平方再减一,即得不均匀度的平方.此处使用不均匀度的平方是因为可以消去绝对值运算和开方运算,减少运算次数、同时也使目标函数变得可导.三、线圈分布规律猜测根据之前的推导易于得出,对一个中心位于柱坐标系中z=h处、半径为R、通电流为I的线圈,对于轴线上、高度为z处的场点,磁感应强度z向分量为由于该函数是关于z的偶函数,奇数阶导数一定都恒等于0,对其求偶数阶导数如下.而对于4个线圈,分别按照题目中定义的参数如下表:线圈1234中心位置-h2-h1h1h2半径R2R1R1R2电流I2I1I1I2可以求得4线圈同时作用时,原点处磁感应强度的z向分量各阶导数为<1><2>而对同一个上述线圈,仍根据之前得出的式子,可以推导出对于与轴线相距x、高度为0处的场点,磁感应强度z向分量为由于问题本身的物理背景,仍然可以看出此函数是关于x的偶函数,奇数阶导数一定都恒等于0.并且根据高等数学的知识,可以证明此函数对于x的各阶导数相当于积分内的函数先对x求偏导,再做对θ的积分.求导工作量很大,借助MATLAB的符号工具箱对其求原点<0,0>处偶数阶导数的值,如下所示.仍与刚才一样代入4个线圈的参数,得到4个线圈共同作用时原点处的各阶导数如下<3><4>至此,得到了<1>~<4>式.令它们均等于0,可以联立求解出4个未知变量R2、h1、h2、I2.联立的方程组如下.得到一组解：h1 = 0.293332954474377I2 =3.000000000000000可以认为按照此规律分布的磁场比较均匀,故可以用此作为优化初值.四、仿真程序编写1.圆周上小段直电流产生的磁场函数dBz<θ,R,r>在柱坐标系中,按照上述推导,编写了在处长为Rdθ的一段电流,方向为,对于位于处的场点,产生的磁场大小在z轴上分量的计算函数如下〔没有乘系数μ/4π〕.dBz = >./sqrt<+R.^2+z.^2>.^3>2.一个圆产生的磁感应强度函数Bz_single在柱坐标系中对位于处的场点,在原点处的半径为R的线圈产生的磁场z轴分量计算函数如下.由于每次计算磁场都会乘系数μ/4π,这部分可以留到函数外面计算,故函数的返回值实际上是没有乘系数μ/4π的磁场强度.函数采用自己编写的梯形公式进行积分,其中step为积分步长,last为上一积分步的被积函数值,now为当前积分步的被积函数值.functionBz = Bz_single% to calculate Bz*4*pi/mu of point in magnetic field,% with integral steps vector thetadBz = >./sqrt<+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; % to store dBz of this stepBz=Bz+; last=now;Bz=I*R*0.5*step*Bz;end对此函数进行正确性测试.在《大学物理学·电磁学》中有一道习题17.6,求一个半圆形电流I在半圆轴线上距离圆心z处的磁场强度.习题答案中给出的磁场强度轴线方向分量为,为了适应此函数,也即.由于题目中是半圆,相当于积分只从0积分到π,在不同的step下得到的积分结果如下表stepπ/1000π/100π/10准确值积分结果z=03.1415926535897933.1415926535897943.1415926535897933.141592653589793z=11.1107207345396001.1107207345395911.1107207345395921.110720734539592z=20.2809925892416300.2809925892416300.2809925892416290.280992589241629z=30.0993458826579590.0993458826579610.0993458826579610.099345882657961可以看出首先函数编写的是准确的,但在这种情况下,积分步长对函数值变化的影响并不是很大,甚至出现了积分步长越粗大,结果反而更精确的反常现象.这有可能是因为π除以一个大数之后变得不精确的原因.为了找到一个经济合适的积分步长,之后再取一个更有一般性的例子,如线圈半径R=1,电流I=2,场点距轴线距离r=0.3,场点距圆心高度z=0.5,这个例子中的各个参数和最终要解决的问题是一个数量级的.测试不同步长情况下的原函数〔从0积分到2π〕计算一个线圈磁感应强度的积分结果如下表stepπ/10000π/1000π/100π/20π/10π/5积分结果8.9576615908211028.9576615908211028.9576615908210998.9576615908210998.9576615908311958.957675492861105误差<1e-15<1e-151e-141e-141e-111e-4此处对θ的步长选取相当于是将圆等效成了多边形,可以看出π/20、π/10的步长已经很能满足一般的计算精度要求了,也很显然可以从直观上感觉,40边形和20边形对于圆的近似程度是很高的.因此最后的程序中采用步长π/20进行计算.再对函数进行效率优化.由于此积分函数在最终的程序中会被反复调用,此函数耗时一定要短.最初想用MATLAB自带的高斯-拉布拉托积分quad<>或quadl<>,但由于其中涉与高阶差分和自适应步长等过程,精度虽高,却不如简单的梯形积分快.且由之前过程可看出,只要步长合适,梯形积分也可以满足精度要求.此函数中在自己编写梯形积分时,每个θ对应的被积函数仅计算一遍,相同的乘法部分均提出到for循环外,没有任何冗余运算,同时也在保证效率前提下将占用变量〔即系统内存〕尽可能的进行了缩减.3.两对圆电流总磁感应强度函数Bz_total其中R2为原图中外侧线圈半径,h1、h2分别为内、外侧线圈到原点的距离,I2为外侧线圈电流,以和I1同向为正,r、z分别为采样点的半径坐标向量、纵坐标向量.本函数4次调用单个圆电流产生磁场的函数,在柱坐标系中将4个线圈的总效应相加,调用一次即可计算所有采样点的z向磁感应强度,结果返回给矩阵变量Bz.编写的程序如下.functionBz = Bz_totalR1 = 1; I1 = 1;% if use parallel computing[X,Z]=meshgrid;Z1 = Z+h2; % coil 1Z2 = Z+h1; % coil 2Z3 = Z-h1; % coil 3Z4 = Z-h2; % coil 4Bz1 = Bz_single;Bz2 = Bz_single;Bz3 = Bz_single;Bz4 = Bz_single; Bz = Bz1+Bz2+Bz3+Bz4;end 值得一提的是,本函数没有采用大家普遍采用的for循环方式来计算400个采样点,而是采用了meshgrid函数生成400个点的坐标矩阵网格,然后对坐标矩阵进行操作,并行的同时计算400个采样点.在不利用其他并行工具的前提下,MATLAB中for循环是一种低效率的串行方式,而通过上网查找资料发现较高版本的MATLAB中矩阵三角函数运算、矩阵开方运算、矩阵点乘除乘方运算都是自动将运算量分配到多个计算单元上并行执行的,能够显著提速.由于我们用到的运算也就是以上几种的组合,因此本函数中先生成了关于r、z坐标向量的两个二维网格矩阵,再直接对这两个网格矩阵进行运算,以达到并行计算的目的.经过实测,其他条件相同的情况下,采用for循环方式做一次优化至3分钟后仍未出结果〔其他组同学的程序平均时间一般也在10分钟左右〕；采用坐标网格矩阵的并行方式做一次优化只需要2秒多.4.不均匀度函数differ由于题中所给采样区的磁场分布一定是对称的,故此处不均匀度函数仅选取了题目中一半仿真区域11×21个点,半径方向20等分取21点,纵轴方向10等分取11点.函数返回值d为所有采样点不均匀度的平方〔这对于找其最小值并无影响〕,输入变量x为一个4维列向量,分别是4个待定参数R2、h1、h2、I2.按照之前推导的简化计算量的公式,编写程序如下.function d = differ% 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*h1:0.8*h1; % z向仅对上半部分进行采样Bz_half = Bz_total;Bz = [Bz_half;Bz_half<2:end,:>];d = mean2/mean2^2-1;end因为此函数中没有绝对值,之前的其他函数部分也没有超越函数,所以至此,待优化的目标函数是可导的,有利于后续优化算法的选择.5.优化主函数选用优化工具箱中的函数fmincon进行优化。
中级消防设施操作员监控26道线下抽考题.pdf
