Helmholtz線圈磁場(chǎng)的仿真研究

1、一、仿真題目要求為了獲得一定區(qū)域上的勻強(qiáng)磁場(chǎng)，可采用多組Helmholtz線圈結(jié)構(gòu)。一種兩對(duì)線圈的結(jié)構(gòu)如下列圖所示。線圈半徑a1，a2，線圈間距h1，h2，以及線圈中通過(guò)電流i1，i2可變化量，如圖1a所示。為了定量衡量關(guān)注區(qū)域的磁場(chǎng)均壓程度，過(guò)軸線做截面ABo1o2，取CD=0.8×AB和Eo3 = 0.8×o1o2，在CD和Eo3線段上每邊均勻取20采樣點(diǎn)，從而形成如圖1b所示的采樣節(jié)點(diǎn)，定義z方向B的不均壓系數(shù)為：其中，為所有采樣點(diǎn)的z方向磁感應(yīng)強(qiáng)度平均值；為第n個(gè)采樣點(diǎn)的z方向磁感應(yīng)強(qiáng)度值，N為采樣點(diǎn)總數(shù)。如給定a1 = 1 m， i1=1 A，且h1 > 0

2、.2 m， 情況下，如果要使得到達(dá)最小，線圈半徑a2，線圈間距h1，h2，以及線圈中通過(guò)電流i2應(yīng)如何取值？(a)線圈結(jié)構(gòu)示意圖b磁場(chǎng)采樣節(jié)點(diǎn)示意圖二、公式推導(dǎo)1.對(duì)于單線圈情況推導(dǎo)磁感應(yīng)強(qiáng)度由于仿真要求在xoz平面上取點(diǎn)，對(duì)磁感應(yīng)強(qiáng)度z軸分量進(jìn)行不均壓度分析，那么只需關(guān)注在xoz平面上點(diǎn)的，推導(dǎo)過(guò)程如下：如上圖所示，位于處的電流元，對(duì)場(chǎng)點(diǎn)P(r,z)處產(chǎn)生的磁感應(yīng)強(qiáng)度z軸分量那么整個(gè)線圈電流元中從0到，進(jìn)行積分得到整個(gè)線圈對(duì)場(chǎng)點(diǎn)P產(chǎn)生的磁感應(yīng)強(qiáng)度z軸分量2.不均勻度公式的簡(jiǎn)化題目中給定的不均勻度為所有采樣點(diǎn)的磁感應(yīng)強(qiáng)度標(biāo)準(zhǔn)差除以平均值的絕對(duì)值，而由數(shù)學(xué)推導(dǎo)可得如下化簡(jiǎn)形式，即采樣點(diǎn)磁感應(yīng)強(qiáng)度

3、平方的均值除以均值的平方再減一，即得不均勻度的平方。此處使用不均勻度的平方是因?yàn)榭梢韵ソ^對(duì)值運(yùn)算和開(kāi)方運(yùn)算，減少運(yùn)算次數(shù)、同時(shí)也使目標(biāo)函數(shù)變得可導(dǎo)。三、線圈分布規(guī)律猜想根據(jù)之前的推導(dǎo)易于得出，對(duì)一個(gè)中心位于柱坐標(biāo)系中z=h處、半徑為R、通電流為I的線圈，對(duì)于軸線上、高度為z處的場(chǎng)點(diǎn)，磁感應(yīng)強(qiáng)度z向分量為由于該函數(shù)是關(guān)于z的偶函數(shù)，奇數(shù)階導(dǎo)數(shù)一定都恒等于0，對(duì)其求偶數(shù)階導(dǎo)數(shù)如下。而對(duì)于4個(gè)線圈，分別按照題目中定義的參數(shù)如下表:線圈1234中心位置-h2-h1h1h2半徑R2R1R1R2電流I2I1I1I2可以求得4線圈同時(shí)作用時(shí)，原點(diǎn)處磁感應(yīng)強(qiáng)度的z向分量各階導(dǎo)數(shù)為 而對(duì)同一個(gè)上述線圈，仍根據(jù)

4、之前得出的式子，可以推導(dǎo)出對(duì)于與軸線相距x、高度為0處的場(chǎng)點(diǎn)，磁感應(yīng)強(qiáng)度z向分量為由于問(wèn)題本身的物理背景，仍然可以看出此函數(shù)是關(guān)于x的偶函數(shù)，奇數(shù)階導(dǎo)數(shù)一定都恒等于0。并且根據(jù)高等數(shù)學(xué)的知識(shí)，可以證明此函數(shù)對(duì)于x的各階導(dǎo)數(shù)相當(dāng)于積分內(nèi)的函數(shù)先對(duì)x求偏導(dǎo)，再做對(duì)的積分。求導(dǎo)工作量很大，借助MATLAB的符號(hào)工具箱對(duì)其求原點(diǎn)(0,0)處偶數(shù)階導(dǎo)數(shù)的值，如下所示。仍與剛剛一樣代入4個(gè)線圈的參數(shù)，得到4個(gè)線圈共同作用時(shí)原點(diǎn)處的各階導(dǎo)數(shù)如下 至此，得到了(1)(4)式。令它們均等于0，可以聯(lián)立求解出4個(gè)未知變量R2、h1、h2、I2。聯(lián)立的方程組如下。得到一組解：R2 = 2.1535923785202

5、54h1 = 0.293332954474377h2 = 1.967170747875403I2 = 3.000000000000000可以認(rèn)為按照此規(guī)律分布的磁場(chǎng)比擬均勻，故可以用此作為優(yōu)化初值。四、仿真程序編寫1.圓周上小段直電流產(chǎn)生的磁場(chǎng)函數(shù)dBz(,R,r)在柱坐標(biāo)系中，按照上述推導(dǎo)，編寫了在(R,)處長(zhǎng)為Rd的一段電流，方向?yàn)?，?duì)于位于(r,z)處的場(chǎng)點(diǎn)，產(chǎn)生的磁場(chǎng)大小在z軸上分量的計(jì)算函數(shù)如下沒(méi)有乘系數(shù)/4。dBz = (theta) (R-r.*cos(theta)./sqrt(r.2-2.*r.*R.*cos(theta)+R.2+z.2).3)2.一個(gè)圓產(chǎn)生的磁感應(yīng)強(qiáng)度函數(shù)B

6、z_single(r,z,R,I)在柱坐標(biāo)系中對(duì)位于(r,z)處的場(chǎng)點(diǎn)，在原點(diǎn)處的半徑為R的線圈產(chǎn)生的磁場(chǎng)z軸分量計(jì)算函數(shù)如下。由于每次計(jì)算磁場(chǎng)都會(huì)乘系數(shù)/4，這局部可以留到函數(shù)外面計(jì)算，故函數(shù)的返回值實(shí)際上是沒(méi)有乘系數(shù)/4的磁場(chǎng)強(qiáng)度。函數(shù)采用自己編寫的梯形公式進(jìn)行積分，其中step為積分步長(zhǎng)，last為上一積分步的被積函數(shù)值，now為當(dāng)前積分步的被積函數(shù)值。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對(duì)此函數(shù)進(jìn)行正確性測(cè)試。在?大學(xué)物理學(xué)(A版)·電磁學(xué)?中有一道習(xí)題1

8、7.6，求一個(gè)半圓形電流I在半圓軸線上距離圓心z處的磁場(chǎng)強(qiáng)度。習(xí)題答案中給出的磁場(chǎng)強(qiáng)度軸線方向分量為，為了適應(yīng)此函數(shù)，也即。由于題目中是半圓，相當(dāng)于積分只從0積分到，在不同的step下得到的積分結(jié)果如下表step/1000/100/10準(zhǔn)確值積分結(jié)果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可以看出首先函數(shù)編寫的是準(zhǔn)確的，但在這種情況下，積分步長(zhǎng)對(duì)函數(shù)值變化的影響并不是很大，甚至出現(xiàn)了積分步長(zhǎng)越粗大，結(jié)果反而更精確的反?，F(xiàn)象。這有可能是因?yàn)槌砸粋€(gè)大數(shù)之后變得不精確的原因。為了找到一個(gè)經(jīng)濟(jì)適宜的積分步長(zhǎng)，之后再取一個(gè)更有一般性的例子，如線圈半徑R=1，電流I=2，場(chǎng)點(diǎn)距軸線距離r=0.3，場(chǎng)點(diǎn)距圓心高度

10、z=0.5，這個(gè)例子中的各個(gè)參數(shù)和最終要解決的問(wèn)題是一個(gè)數(shù)量級(jí)的。測(cè)試不同步長(zhǎng)情況下的原函數(shù)從0積分到2計(jì)算一個(gè)線圈磁感應(yīng)強(qiáng)度的積分結(jié)果如下表step/10000/1000/100/20/10/5積分結(jié)果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此處對(duì)的步長(zhǎng)選取相當(dāng)于是將圓等效成了多邊形，可以看出/20、/10的步長(zhǎng)已經(jīng)很能滿足一般的

11、計(jì)算精度要求了，也很顯然可以從直觀上感覺(jué)，40邊形和20邊形對(duì)于圓的近似程度是很高的。因此最后的程序中采用步長(zhǎng)/20進(jìn)行計(jì)算。再對(duì)函數(shù)進(jìn)行效率優(yōu)化。由于此積分函數(shù)在最終的程序中會(huì)被反復(fù)調(diào)用，此函數(shù)耗時(shí)一定要短。最初想用MATLAB自帶的高斯-拉布拉托積分quad()或quadl()，但由于其中涉及高階差分和自適應(yīng)步長(zhǎng)等過(guò)程，精度雖高，卻不如簡(jiǎn)單的梯形積分快。且由之前過(guò)程可看出，只要步長(zhǎng)適宜，梯形積分也可以滿足精度要求。此函數(shù)中在自己編寫梯形積分時(shí)，每個(gè)對(duì)應(yīng)的被積函數(shù)僅計(jì)算一遍，相同的乘法局部均提出到for循環(huán)外，沒(méi)有任何冗余運(yùn)算，同時(shí)也在保證效率前提下將占用變量即系統(tǒng)內(nèi)存盡可能的進(jìn)行了縮減。3

12、.兩對(duì)圓電流總磁感應(yīng)強(qiáng)度函數(shù)Bz_total(R2,h1,h2,I2,r,z)其中R2為原圖中外側(cè)線圈半徑，h1、h2分別為內(nèi)、外側(cè)線圈到原點(diǎn)的距離，I2為外側(cè)線圈電流，以和I1同向?yàn)檎?，r、z分別為采樣點(diǎn)的半徑坐標(biāo)向量、縱坐標(biāo)向量。本函數(shù)4次調(diào)用單個(gè)圓電流產(chǎn)生磁場(chǎng)的函數(shù)，在柱坐標(biāo)系中將4個(gè)線圈的總效應(yīng)相加，調(diào)用一次即可計(jì)算所有采樣點(diǎn)的z向磁感應(yīng)強(qiáng)度，結(jié)果返回給矩陣變量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 值得一提的是，本函數(shù)沒(méi)有采用大家普遍采用的for循環(huán)方式來(lái)計(jì)算400個(gè)采樣點(diǎn)，而是采用了meshgrid函數(shù)生成400個(gè)點(diǎn)的坐標(biāo)矩陣網(wǎng)格，然后對(duì)坐標(biāo)矩陣進(jìn)行操作，并行的同時(shí)計(jì)算40

14、0個(gè)采樣點(diǎn)。在不利用其他并行工具的前提下，MATLAB中for循環(huán)是一種低效率的串行方式，而通過(guò)上網(wǎng)查找資料發(fā)現(xiàn)較高版本的MATLAB中矩陣三角函數(shù)運(yùn)算、矩陣開(kāi)方運(yùn)算、矩陣點(diǎn)乘除乘方運(yùn)算都是自動(dòng)將運(yùn)算量分配到多個(gè)計(jì)算單元上并行執(zhí)行的，能夠顯著提速。由于我們用到的運(yùn)算也就是以上幾種的組合，因此本函數(shù)中先生成了關(guān)于r、z坐標(biāo)向量的兩個(gè)二維網(wǎng)格矩陣，再直接對(duì)這兩個(gè)網(wǎng)格矩陣進(jìn)行運(yùn)算，以到達(dá)并行計(jì)算的目的。經(jīng)過(guò)實(shí)測(cè)，其他條件相同的情況下，采用for循環(huán)方式做一次優(yōu)化至3分鐘后仍未出結(jié)果其他組同學(xué)的程序平均時(shí)間一般也在10分鐘左右；采用坐標(biāo)網(wǎng)格矩陣的并行方式做一次優(yōu)化只需要2秒多。4.不均勻度函數(shù)diff

15、er(x)由于題中所給采樣區(qū)的磁場(chǎng)分布一定是對(duì)稱的，故此處不均勻度函數(shù)僅選取了題目中一半仿真區(qū)域11×21個(gè)點(diǎn)，半徑方向20等分取21點(diǎn)，縱軸方向10等分取11點(diǎn)。函數(shù)返回值d為所有采樣點(diǎn)不均勻度的平方這對(duì)于找其最小值并無(wú)影響，輸入變量x為一個(gè)4維列向量，分別是4個(gè)待定參數(shù)R2、h1、h2、I2。按照之前推導(dǎo)的簡(jiǎn)化計(jì)算量的公式，編寫程序如下。function d = differ(x)% x=R2;h1;h2;I2% d 為不均勻系數(shù)delta的平方R2=x(1);h1=x(2);h2=x(3);I2=x(4);% 每邊采樣點(diǎn)數(shù)設(shè)置r = 0:0.04:0.8;z = 0:0.08*

16、h1:0.8*h1; % z向僅對(duì)上半局部進(jìn)行采樣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因?yàn)榇撕瘮?shù)中沒(méi)有絕對(duì)值，之前的其他函數(shù)局部也沒(méi)有超越函數(shù)，所以至此，待優(yōu)化的目標(biāo)函數(shù)是可導(dǎo)的，有利于后續(xù)優(yōu)化算法的選擇。5.優(yōu)化主函數(shù)選用優(yōu)化工具箱中的函數(shù)fmincon進(jìn)行優(yōu)化，找到目標(biāo)函數(shù)的最小值。其中約束條件的設(shè)置在矩陣A1中和向量v1、v2中，限定了h2>h1、h2>0.2、-3<I2<3。向量v1、v2的其他局部為人為

17、設(shè)定的另外約束，為了方便找到最優(yōu)解。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; 五、數(shù)值試驗(yàn)1.按之前導(dǎo)數(shù)為零的推導(dǎo)選取優(yōu)化初值首先選取了之前按照原點(diǎn)處各階導(dǎo)數(shù)為0推導(dǎo)出的一組參數(shù)作為初值，變量范圍設(shè)在其附近，得到了一組優(yōu)化結(jié)果，畫出其磁感應(yīng)強(qiáng)度分布如下列圖(a)。但隨后發(fā)現(xiàn)當(dāng)初值取為其他值時(shí)會(huì)出現(xiàn)不均勻度更小的解，如下列圖(b)。這說(shuō)明之前的推導(dǎo)是不正確的，它僅僅保證了原點(diǎn)附近磁場(chǎng)足夠均勻，卻不能保證整個(gè)采樣區(qū)域的不均勻度最小。(a)按照導(dǎo)數(shù)為零得到的結(jié)果(b)按照另一組初值得到的結(jié)果2.按數(shù)值試驗(yàn)逐漸縮小范圍的方法選取初值我們希望能通過(guò)數(shù)值實(shí)驗(yàn)先計(jì)算幾組參數(shù)值下的不均勻系數(shù)平方，再進(jìn)一步縮小范圍尋找其最小值，確定所優(yōu)化目標(biāo)的大致位置。由于要優(yōu)化四個(gè)參數(shù)，假設(shè)參數(shù)值一起改變，不太方便觀察，可先固定兩個(gè)參數(shù)不變，改變其他兩個(gè)，再反過(guò)來(lái)改變另兩個(gè)即可。根據(jù)之前幾次優(yōu)化發(fā)現(xiàn)，一般最后結(jié)果都會(huì)有h1=0.2，I2=3，而另外兩個(gè)量