Helmholtz線圈磁場的仿真研究模板

資料內(nèi)容僅供您學習參考，如有不當或者侵權，請聯(lián)系改正或者刪除。一、仿真題目要求為了獲得一定區(qū)域上的勻強磁場,可采用多組Helmholtz線圈結構。一種兩對線圈的結構如下圖所示。線圈半徑a1,a2,線圈間距h1,h2,以及線圈中經(jīng)過電流i1,i2可變化量,如圖1(a)所示。為了定量衡量關注區(qū)域的磁場均壓程度,過軸線做截面ABo1o2,取CD=0.8×AB和Eo3=0.8×o1o2,在CD和Eo3線段上每邊均勻取20采樣點,從而形成如圖1(b)所示的采樣節(jié)點,定義z方向B的不均壓系數(shù)為:其中,為所有采樣點的z方向磁感應強度平均值;為第n個采樣點的z方向磁感應強度值,N為采樣點總數(shù)。如給定a1=1m,i1=1A,且h1>0.2m,情況下,如果要使得δ達到最小,線圈半徑a2,線圈間距h1,h2,以及線圈中經(jīng)過電流i2應如何取值?(a)線圈結構示意圖(b)磁場采樣節(jié)點示意圖二、公式推導1.對于單線圈情況推導磁感應強度由于仿真要求在xoz平面上取點,對磁感應強度z軸分量進行不均壓度分析,則只需關注在xoz平面上點的,推導過程如下:如上圖所示,位于處的電流元,對場點P(r,z)處產(chǎn)生的磁感應強度z軸分量則整個線圈電流元中從0到,進行積分得到整個線圈對場點P產(chǎn)生的磁感應強度z軸分量 2.不均勻度公式的簡化題目中給定的不均勻度為所有采樣點的磁感應強度標準差除以平均值的絕對值,而由數(shù)學推導可得如下化簡形式,即采樣點磁感應強度平方的均值除以均值的平方再減一,即得不均勻度的平方。此處使用不均勻度的平方是因為能夠消去絕對值運算和開方運算,減少運算次數(shù)、同時也使目標函數(shù)變得可導。三、線圈分布規(guī)律猜測根據(jù)之前的推導易于得出,對一個中心位于柱坐標系中z=h處、半徑為R、通電流為I的線圈,對于軸線上、高度為z處的場點,磁感應強度z向分量為由于該函數(shù)是關于z的偶函數(shù),奇數(shù)階導數(shù)一定都恒等于0,對其求偶數(shù)階導數(shù)如下。而對于4個線圈,分別按照題目中定義的參數(shù)如下表:線圈1234中心位置-h2-h1h1h2半徑R2R1R1R2電流I2I1I1I2能夠求得4線圈同時作用時,原點處磁感應強度的z向分量各階導數(shù)為而對同一個上述線圈,仍根據(jù)之前得出的式子,能夠推導出對于與軸線相距x、高度為0處的場點,磁感應強度z向分量為由于問題本身的物理背景,依然能夠看出此函數(shù)是關于x的偶函數(shù),奇數(shù)階導數(shù)一定都恒等于0。而且根據(jù)高等數(shù)學的知識,能夠證明此函數(shù)對于x的各階導數(shù)相當于積分內(nèi)的函數(shù)先對x求偏導,再做對θ的積分。求導工作量很大,借助MATLAB的符號工具箱對其求原點(0,0)處偶數(shù)階導數(shù)的值,如下所示。仍與剛才一樣代入4個線圈的參數(shù),得到4個線圈共同作用時原點處的各階導數(shù)如下至此,得到了(1)~(4)式。令它們均等于0,能夠聯(lián)立求解出4個未知變量R2、h1、h2、I2。聯(lián)立的方程組如下。得到一組解:R2=2.0254h1=0.377h2=1.9673I2=3.000能夠認為按照此規(guī)律分布的磁場比較均勻,故能夠用此作為優(yōu)化初值。四、仿真程序編寫1.圓周上小段直電流產(chǎn)生的磁場函數(shù)dBz(θ,R,r)在柱坐標系中,按照上述推導,編寫了在(R,θ)處長為Rdθ的一段電流,方向為,對于位于(r,z)處的場點,產(chǎn)生的磁場大小在z軸上分量的計算函數(shù)如下(沒有乘系數(shù)μ/4π)。dBz=@(theta)(R-r.*cos(theta))./sqrt((r.^2-2.*r.*R.*cos(theta)+R.^2+z.^2).^3)2.一個圓產(chǎn)生的磁感應強度函數(shù)Bz_single(r,z,R,I)在柱坐標系中對位于(r,z)處的場點,在原點處的半徑為R的線圈產(chǎn)生的磁場z軸分量計算函數(shù)如下。由于每次計算磁場都會乘系數(shù)μ/4π,這部分能夠留到函數(shù)外面計算,故函數(shù)的返回值實際上是沒有乘系數(shù)μ/4π的磁場強度。函數(shù)采用自己編寫的梯形公式進行積分,其中step為積分步長,last為上一積分步的被積函數(shù)值,now為當前積分步的被積函數(shù)值。functionBz=Bz_single(r,z,R,I)%tocalculateBz*4*pi/muofpoint(x,z)inmagneticfield,%withintegralstepsvectorthetadBz=@(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);%tostoredBzoflaststepfortheta=step:step:2*pi;now=dBz(theta);%tostoredBzofthisstepBz=Bz+(now+last);last=now;Bz=I*R*0.5*step*Bz;end對此函數(shù)進行正確性測試。在《大學物理學(A版)·電磁學》中有一道習題17.6,求一個半圓形電流I在半圓軸線上距離圓心z處的磁場強度。習題答案中給出的磁場強度軸線方向分量為,為了適應此函數(shù),也即。由于題目中是半圓,相當于積分只從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能夠看出首先函數(shù)編寫的是準確的,但在這種情況下,積分步長對函數(shù)值變化的影響并不是很大,甚至出現(xiàn)了積分步長越粗大,結果反而更精確的反?，F(xiàn)象。這有可能是因為π除以一個大數(shù)之后變得不精確的原因。為了找到一個經(jīng)濟合適的積分步長,之后再取一個更有一般性的例子,如線圈半徑R=1,電流I=2,場點距軸線距離r=0.3,場點距圓心高度z=0.5,這個例子中的各個參數(shù)和最終要解決的問題是一個數(shù)量級的。測試不同步長情況下的原函數(shù)(從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的步長已經(jīng)很能滿足一般的計算精度要求了,也很顯然能夠從直觀上感覺,40邊形和20邊形對于圓的近似程度是很高的。因此最后的程序中采用步長π/20進行計算。再對函數(shù)進行效率優(yōu)化。由于此積分函數(shù)在最終的程序中會被重復調(diào)用,此函數(shù)耗時一定要短。最初想用MATLAB自帶的高斯-拉布拉托積分quad()或quadl(),但由于其中涉及高階差分和自適應步長等過程,精度雖高,卻不如簡單的梯形積分快。且由之前過程可看出,只要步長合適,梯形積分也能夠滿足精度要求。此函數(shù)中在自己編寫梯形積分時,每個θ對應的被積函數(shù)僅計算一遍,相同的乘法部分均提出到for循環(huán)外,沒有任何冗余運算,同時也在保證效率前提下將占用變量(即系統(tǒng)內(nèi)存)盡可能的進行了縮減。3.兩對圓電流總磁感應強度函數(shù)Bz_total(R2,h1,h2,I2,r,z)其中R2為原圖中外側(cè)線圈半徑,h1、h2分別為內(nèi)、外側(cè)線圈到原點的距離,I2為外側(cè)線圈電流,以和I1同向為正,r、z分別為采樣點的半徑坐標向量、縱坐標向量。本函數(shù)4次調(diào)用單個圓電流產(chǎn)生磁場的函數(shù),在柱坐標系中將4個線圈的總效應相加,調(diào)用一次即可計算所有采樣點的z向磁感應強度,結果返回給矩陣變量Bz。編寫的程序如下。functionBz=Bz_total(R2,h1,h2,I2,r,z)R1=1;I1=1;%ifuseparallelcomputing[X,Z]=meshgrid(r,z);Z1=Z+h2;%coil1Z2=Z+h1;%coil2Z3=Z-h1;%coil3Z4=Z-h2;%coil4Bz1=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ù)沒有采用大家普遍采用的for循環(huán)方式來計算400個采樣點,而是采用了meshgrid函數(shù)生成400個點的坐標矩陣網(wǎng)格,然后對坐標矩陣進行操作,并行的同時計算400個采樣點。在不利用其它并行工具的前提下,MATLAB中for循環(huán)是一種低效率的串行方式,而經(jīng)過上網(wǎng)查找資料發(fā)現(xiàn)較高版本的MATLAB中矩陣三角函數(shù)運算、矩陣開方運算、矩陣點乘除乘方運算都是自動將運算量分配到多個計算單元上并行執(zhí)行的,能夠顯著提速。由于我們用到的運算也就是以上幾種的組合,因此本函數(shù)中先生成了關于r、z坐標向量的兩個二維網(wǎng)格矩陣,再直接對這兩個網(wǎng)格矩陣進行運算,以達到并行計算的目的。經(jīng)過實測,其它條件相同的情況下,采用for循環(huán)方式做一次優(yōu)化至3分鐘后仍未出結果(其它組同學的程序平均時間一般也在10分鐘左右);采用坐標網(wǎng)格矩陣的并行方式做一次優(yōu)化只需要2秒多。4.不均勻度函數(shù)differ(x)由于題中所給采樣區(qū)的磁場分布一定是對稱的,故此處不均勻度函數(shù)僅選取了題目中一半仿真區(qū)域11×21個點,半徑方向20等分取21點,縱軸方向10等分取11點。函數(shù)返回值d為所有采樣點不均勻度的平方(這對于找其最小值并無影響),輸入變量x為一個4維列向量,分別是4個待定參數(shù)R2、h1、h2、I2。按照之前推導的簡化計算量的公式,編寫程序如下。functiond=differ(x)%x=[R2;h1;h2;I2]%d為不均勻系數(shù)delta的平方R2=x(1);h1=x(2);h2=x(3);I2=x(4);%每邊采樣點數(shù)設置r=0:0.04:0.8;z=0:0.08*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因為此函數(shù)中沒有絕對值,之前的其它函數(shù)部分也沒有超越函數(shù),因此至此,待優(yōu)化的目標函數(shù)是可導的,有利于后續(xù)優(yōu)化算法的選擇。5.優(yōu)化主函數(shù)選用優(yōu)化工具箱中的函數(shù)fmincon進行優(yōu)化,找到目標函數(shù)的最小值。其中約束條件的設置在矩陣A1中和向量v1、v2中,限定了h2>h1、h2>0.2、-3<I2<3。向量v1、v2的其它部分為人為設定的另外約束,為了方便找到最優(yōu)解。Clearall;clc;tic;mu=4*pi*1e-7;R1=1;I1=1;%variablex=[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)toc;五、數(shù)值試驗1.按之前導數(shù)為零的推導選取優(yōu)化初值首先選取了之前按照原點處各階導數(shù)為0推導出的一組參數(shù)作為初值,變量范圍設在其附近,得到了一組優(yōu)化結果,畫出其磁感應強度分布如下圖(a)。但隨后發(fā)現(xiàn)當初值取為其它值時會出現(xiàn)不均勻度更小的解,如下圖(b)。這說明之前的推導是不正確的,它僅僅保證了原點附近磁場足夠均勻,卻不能保證整個采樣區(qū)域的不均勻度最小。(a)按照導數(shù)為零得到的結果(b)按照另一組初值得到的結果2.按數(shù)值試驗逐漸縮小范圍的方法選取初值我們希望能經(jīng)過數(shù)值實驗先計算幾組參數(shù)值下的不均勻系數(shù)平方,再進一步縮小范圍尋找其最小值,確定所優(yōu)化目標的大致位置。由于要優(yōu)化四個參數(shù),若參數(shù)值一起改變,不太方便觀察,可先固定兩個參數(shù)不變,改變其它兩個,再反過來改變另兩個即可。根據(jù)之前幾次優(yōu)