基于矩量法的二維金屬體散射(內(nèi)含matlab程序)

1、基于矩量法的二維金屬體散射計(jì)算1 問(wèn)題的描述本題是用矩量法計(jì)算二維金屬圓柱體的散射場(chǎng)，如圖所示為一圓柱體和一個(gè)橢圓柱的截面，為了計(jì)算簡(jiǎn)單，選入射波為垂直z軸入射的TM或TE平面波 y 2x 2 矩量法求解過(guò)程2.1 電場(chǎng)積分方程 2.1.1問(wèn)題的分析由麥克斯韋方程組 (1) (2)可得電場(chǎng)積分方程為 (3)表示在圓柱表面的面電流在遠(yuǎn)處產(chǎn)生的總場(chǎng)。設(shè)入射場(chǎng)為E，散射場(chǎng)為E，由金屬表面的邊界條件=0 (4)得 (5)2.1.2 離散化設(shè)入射波為,將散射體截面C分為N份C，用點(diǎn)匹配法對(duì)上述積分式子進(jìn)行離散化， 即基函數(shù)可取 (6)可得下列離散方程：PJ=b (7)其中： (8) (9)當(dāng)mn時(shí)， (

2、10)當(dāng)m=n時(shí) 解析積分為 (11)其中=1.781,e=2.7182.1.3方程組的求解可用LU分解求解方程組，即P=LU ，其中P為可逆矩陣，L為上三角矩陣，U為下三角矩陣，則可利用這兩個(gè)基本的三角矩陣進(jìn)行求解J，求出J之后，就可求散射場(chǎng) (12) (13)與二維場(chǎng)中的散射截面 (14)2.1.4輸出結(jié)果的驗(yàn)證此散射問(wèn)題也可用模式展開(kāi)法進(jìn)行求解，可用此結(jié)果對(duì)本問(wèn)題進(jìn)行驗(yàn)證。所得J為 (15)2.2 磁場(chǎng)積分方程對(duì)于TE波垂直與z方向入射時(shí)的金屬體的散射。對(duì)于一般的TE波而言只有場(chǎng)分量，電流密度方程只有橫向分量。則MFIE為： (16)其中(17) yytnxx (18) 其中t表示邊界上

3、的一點(diǎn)，是和X的夾角。根根據(jù)前面的過(guò)程，圓柱邊界分成N分。等效電流密度可以近似為一些脈沖函數(shù)的迭加： (19)其中 (20)則得到 (21)矩陣非對(duì)角元 (22) (23) 在上認(rèn)為是常量故 (24)對(duì)角元 (25)其中 回波寬度的近似公式為： (26)3 計(jì)算機(jī)數(shù)值實(shí)驗(yàn)及分析本論文通過(guò)數(shù)值計(jì)算驗(yàn)證前面理論分析的結(jié)果，并對(duì)數(shù)值計(jì)算結(jié)果進(jìn)行分析。分別以金屬圓柱體和金屬橢圓體為計(jì)算例子，做數(shù)值實(shí)驗(yàn)和分析。所使用的計(jì)算機(jī)程序是商業(yè)軟件MATLAB6.5，數(shù)值實(shí)驗(yàn)在本人機(jī)子（celeron4 1.8G CPU 128M內(nèi)存）,操作系統(tǒng)是windows xp。3.1 二維金屬圓柱體的散射基于上面的分析，

4、考慮垂直z方向入射的橫向磁波（TM）,離散方程為（7）,編程的基本思路是對(duì)(10)式和（11）式編程實(shí)現(xiàn)，得出p矩陣，再由(9)式得出b列，用MATLAB6.5軟件上的線(xiàn)性方程組直接求解法求解出J。散射截面（回波寬度）可以通過(guò)（14）式離散計(jì)算出來(lái)。計(jì)算例子是一個(gè)z方向均勻且無(wú)限長(zhǎng)的金屬圓柱，半徑為1.5米（），金屬圓柱中心和z軸重合，入射波為z方向極化，幅值為1，從負(fù)x軸方向垂直z軸入射的TM平面波，工作頻率為100MHz，波長(zhǎng)米。由于是金屬體且z方向均勻，可以只考慮對(duì)垂直z軸截面的圓周進(jìn)行剖分并計(jì)算。下圖給出了720個(gè)剖分下電流密度分布的計(jì)算結(jié)果與近似解析解的比較，其中近似解析解是根據(jù)導(dǎo)波

5、理論書(shū)上3.48式(本文的（15）式)在n=-36到n=36下計(jì)算出來(lái)的。計(jì)算時(shí)入射角取為0度。其中x軸為，y軸為電流密度，由圖可見(jiàn)，電流密度分布和近似解析解無(wú)論幅度相位之間都有著非常好的吻合。計(jì)算所得的總等效電流Iz-0.0079 + 0.0083i,而在剖分精度為180時(shí)，計(jì)算所得的總等效電流Iz =-0.0084 + 0.0083i。而解析解的總電流Iz =-0.0077 + 0.0083i，可見(jiàn)隨著剖分精度的增加，計(jì)算結(jié)果收斂于解析解。 圖1（a）EFIE ，剖分精度720 圖1（b）EFIE 近似解析解下圖給出的是回波寬度的分布 圖1（c）EFIE 剖分精度720其中x軸為，y軸為，

6、據(jù)個(gè)人粗略分析應(yīng)該基本符合事實(shí)。由于沒(méi)能得到回波寬度的解析解，沒(méi)能作進(jìn)一步的分析比較。下面給出入射角為90度，半徑，而其他條件不變的情況下，所得的計(jì)算結(jié)果。 （d）EFIE 剖分精度720由此可見(jiàn)，相對(duì)前面那種情況，入射角變化90度，等效電流密度分布也相應(yīng)有90度的相移，回波寬度的幅度減小了很多，但大體的形狀保持不變。這時(shí)候的總電流Iz =0.0067 + 0.0028i。3.2 TM波入射金屬橢圓柱的散射對(duì)于二維金屬橢圓柱體的散射這種情況，由于圓柱體是橢圓柱的特殊情況，所以解題的基本思路基本一樣，就是對(duì)每個(gè)剖分步長(zhǎng)用數(shù)值積分得到，這樣有利于得到精確的計(jì)算結(jié)果。實(shí)踐的過(guò)程也證明了這一點(diǎn)，當(dāng)每個(gè)

7、剖分步長(zhǎng)用兩個(gè)剖分點(diǎn)的直線(xiàn)距離來(lái)近似的話(huà)，帶來(lái)很大的誤差，而用數(shù)值積分得到的結(jié)果和解析解很好地吻合。計(jì)算例子是一個(gè)z方向均勻且無(wú)限長(zhǎng)的橢圓柱，長(zhǎng)軸，短軸，即金屬圓柱中心和z軸重合，即橢圓方程為。入射波為z方向極化，幅值為1，從負(fù)x軸方向垂直z軸入射的TM平面波，即入射角為0度。工作頻率為100MHz，波長(zhǎng)米。由于是金屬體且z方向均勻，可以只考慮對(duì)垂直z軸截面的橢圓周進(jìn)行剖分并計(jì)算。圖3（a）(b)給出了1000個(gè)剖分和2000個(gè)剖分下的電流密度分布的計(jì)算結(jié)果，圖3(c)給出解析解以作為比較，而圖4(a)(b)給出了上述剖分精度下的回波寬度的計(jì)算結(jié)果，作為比較圖4(c)給出回波寬度的解析解。其中

8、解析解來(lái)自參考書(shū)計(jì)算電磁場(chǎng)的矩量法。為了方便與參考書(shū)中的解析解比較，x和y軸的參量都相應(yīng)做了變化。下圖中的x軸S為角度的歸一化，左端為S=0,右端為S=1。Y軸是 圖3(a) EFIE 剖分精度1000圖3(b) EFIE 剖分精度2000 圖3(c) EFIE 剖分精度2700 圖3(d) EFIE 安得列解由上圖可見(jiàn)計(jì)算結(jié)果和參考書(shū)提供的解析解能夠有很好的吻合，而且隨著剖分加細(xì)，結(jié)果更趨接近于解析解。由于本人機(jī)子配置較低，難以對(duì)更多剖分點(diǎn)的情況進(jìn)行運(yùn)算，但可以預(yù)見(jiàn)隨著剖分點(diǎn)的增加，計(jì)算結(jié)果與解析解更好地吻合。但隨著剖分?jǐn)?shù)N的增大，計(jì)算方法所用的近似不能收斂于解析解，這是因?yàn)闀r(shí)的，在極限時(shí)是

9、不正確的。證明了計(jì)算方法的正確性。也可以看到要是用磁場(chǎng)積分方程（MFIE）可以得到更好的解。下圖是回波寬度（散射截面）的方向圖，其中x軸是角度，y軸是。 圖4（a）EFIE 剖分精度1000圖4(b) EFIE 剖分精度2000圖4（c）EFIE回波寬度安得列解比較上圖可得，計(jì)算結(jié)果和解析解幾乎完全一致，可以注意到即使電流密度分布顯著不同，當(dāng)這兩種情況得到的結(jié)果卻是幾乎完全一致的。這是因?yàn)槭请娏鱆的連續(xù)性泛函，因此J在精確值附近的大小變化是不敏感的。3.3 TE波入射金屬橢圓柱的散射與上面同樣條件，把入射波換為T(mén)E波，理論分析可見(jiàn)前面的2.2部分?；谏厦娴姆治觯紤]垂直z方向入射的橫向電波（

10、TE）,離散方程為（21）,編程的基本思路是對(duì)(24)式編程實(shí)現(xiàn)，得出Z矩陣，再由(17)式得出H列，用MATLAB6.5軟件上的線(xiàn)性方程組直接求解法求解出J。散射截面（回波寬度）可以通過(guò)（26）式離散計(jì)算出來(lái)。圖5給出剖分為720份時(shí)的計(jì)算結(jié)果，并給出相應(yīng)的解析解，以資比較。解析解來(lái)自參考書(shū)計(jì)算電磁場(chǎng)的矩量法。 圖5(a) MFIE 剖分精度720 圖5（b）MFIE安得列解可見(jiàn)，計(jì)算結(jié)果與解析解大體上符合，但還是存在較大的差別，原因估計(jì)是剖分精度不夠，還有數(shù)值計(jì)算P矩陣時(shí)引進(jìn)了近似，由于時(shí)間倉(cāng)促，沒(méi)能在離散化方程時(shí)考慮更好的近似方法，這有待于進(jìn)一步的探討和研究。圖6是TE波入射金屬橢圓柱的

11、回波寬度（散射截面）的方向圖，其中x軸是角度，y軸是。可見(jiàn)計(jì)算結(jié)果和解析解很好地符合，可見(jiàn)雖然電流分布計(jì)算結(jié)果和所給的解析解有較大的誤差，但回波寬度的計(jì)算結(jié)果和解析解確幾乎完全一致，這是因?yàn)槭请娏鱆的連續(xù)性泛函，因此J在精確值附近的大小變化是不敏感的。 圖6(a) MFIE 剖分精度720 圖6（b）MFIE安得列解4 存在問(wèn)題和心得存在的問(wèn)題之一：沒(méi)考慮到內(nèi)諧振問(wèn)題，進(jìn)一步的工作應(yīng)該把電場(chǎng)積分方程（EFIE）和磁場(chǎng)積分方程（MFIE）組成聯(lián)合積分方程（CFIE），來(lái)解決這個(gè)問(wèn)題。存在的問(wèn)題之二：使用的軟件MATLAB6.5雖然功能強(qiáng)大，但是運(yùn)算效率不高，需要占用的內(nèi)存大和運(yùn)行時(shí)間較長(zhǎng)，而比不

12、上用C語(yǔ)言或FORTRAN語(yǔ)言編寫(xiě)的程序效率高。沒(méi)有利用到課本所介紹的快速多極子技術(shù)。存在的問(wèn)題之三：沒(méi)有實(shí)現(xiàn)在TM波入射情形下用磁場(chǎng)積分方程（MFIE）計(jì)算散射場(chǎng)，而據(jù)理論分析，用MFIE應(yīng)該能夠得到更好的條件數(shù)，計(jì)算的結(jié)果也能更好地收斂于解析解。這有待于工作的進(jìn)一步深入。存在問(wèn)題之四：限于作者知識(shí)和經(jīng)驗(yàn)的不足，沒(méi)能有更為深厚的理論認(rèn)識(shí)做指導(dǎo)，對(duì)結(jié)果的分析未免有失偏頗。進(jìn)一步的工作應(yīng)該朝著三維散射，介質(zhì)體散射的方向進(jìn)行。心得：本文凝結(jié)著本組成員的心血，期間經(jīng)歷幾多挫折，幸好一一克服了，要說(shuō)有什么心得的話(huà)，第一要對(duì)電磁理論有深刻的認(rèn)識(shí)和理解，第二要有深厚的數(shù)學(xué)基礎(chǔ)，對(duì)電磁積分方程的性能有深入的

13、理解。第三要熟悉熟練掌握編程語(yǔ)言MATLAB，這是實(shí)現(xiàn)的工具。鳴謝：第一應(yīng)該感謝盛新慶老師的悉心指導(dǎo)，第二感激本組成員的通力合作和不懈的努力。5 參考文獻(xiàn) 1 盛新慶. 計(jì)算電磁學(xué)要論. 北京：科學(xué)出版社，2004 2 Roger F.Harrington. Field Compution by Moment Methods. New York:The Macmillan Company, 1968 3 Andrew F.Peterson,Scott L.Ray and Raj Mittra. Computational Methods for Electromagnetics. New Yo

14、rk:IEEE PRESS,1998 4 張志涌等. 精MATLAB6.5版. 北京：北京航空航天大學(xué)出版社，2003 6 程序附錄由于作者已經(jīng)對(duì)程序做了很好的注釋?zhuān)瑧?yīng)該具備matlab基礎(chǔ)的人一般都能夠看懂，所以不準(zhǔn)備做再多的解釋說(shuō)明。只簡(jiǎn)單附在后面，以供參考：5.1金屬圓柱體散射的程序:function win(NU,L) %金屬圓柱體散射的程序，計(jì)算等效表面電流，回波寬度,NU為波長(zhǎng)，L為半徑Z=377; %特性阻抗K=2*pi/NU; %波數(shù)m=720; fai=0; %入射方向和x軸的夾角R=L*NU; %金屬圓柱半徑h=2*pi*R/m; %剖分步長(zhǎng)P=zeros(m,m); %生

15、成Pmn矩陣框架Q=(pi/m:2*pi/m:2*pi*(m-1/2)/m); %角度細(xì)分，分為m分X=R*cos(Q); %x分量長(zhǎng)度，隨角度Q變化Y=R*sin(Q); %y分量長(zhǎng)度，隨角度Q變化dx=zeros(m,m); %生成矩陣 dy=zeros(m,m);for n=1:m %對(duì)每個(gè)x和y分別計(jì)算(x-xm)2和(y-ym)2dx(n,:)=(X-X(n).2;dy(n,:)=(Y-Y(n).2;endI=eye(m);d=dx+dy+I; %(x-xm)2+(y-ym)2,對(duì)角線(xiàn)元素置為1，可以直接用hankel函數(shù)運(yùn)算x=K*sqrt(d); H=besselh(0,2,x)

16、; %hankel functionP=Z*K*h*H/4; %得到P矩陣Pnn=Z*K*h*(1-i*2*log10(1.781*K*h/(4*2.718)/pi)/4; %對(duì)m=n時(shí)，計(jì)算Pnnfor n=1:m P(n,n)=Pnn; %P矩陣對(duì)角元素賦值，所有對(duì)角元素相同endb=exp(-i*K*(X*cos(fai)+Y*sin(fai); %b矩陣,行，要轉(zhuǎn)化為列M=b.; %b array是復(fù)數(shù)，不能用b取得它的轉(zhuǎn)置，而必須用b.。否則得出來(lái)的電流分布有180度的相移j=PM; %求解方程，得到電流密度分布arrayfigure(1)w2=(1:1:m);plot(w2*360

17、/m,abs(j),.) %畫(huà)圖，電流密度在0360度的分布圖grid ontitle(TM波入射的金屬圓柱表面等效電流密度分布圖)xlabel(degrees)ylabel(current dense)Iz=sum(h*j) %總電流w=(0:pi/180:pi*179/180);for n=1:180 %計(jì)算散射截面（回波寬度）g=exp(i*K*(X*cos(w(n)+fai)+Y*sin(w(n)+fai)*h; G=g.*j.; %j array是復(fù)數(shù)，不能用j取得它的轉(zhuǎn)置，而必須用j.。否則得出來(lái)的值有180度的相移T=abs(sum(G);out(n)=K*Z2*T2/4;end

18、W=(0:1:179);figure(2)plot(W,20*log10(out/NU2),.) %畫(huà)圖，回波寬度，0180度grid ontitle(TM波入射金屬圓柱體的回波寬度分布圖)xlabel(degrees)ylabel(echo width/wavelength2/dB)5.2 金屬圓柱體的表面等效電流密度分布圖近似的解析解程序：function current %3.48式數(shù)值計(jì)算,以便對(duì)前面的計(jì)算結(jié)果作檢驗(yàn)。%計(jì)算金屬圓柱體的表面等效電流密度分布圖,得出近似的解析解。%各個(gè)條件與程序win完全一至R=1.5; %金屬圓柱半徑 Z=377; %特性阻抗 K=2*pi/3; %波

19、數(shù)h=2*pi*R/360; %剖分步長(zhǎng)w=(0:2*pi/360:2*pi*359/360); %角度細(xì)分為360分for m=1:360 %對(duì)每個(gè)角度分別計(jì)算Jzfor n=1:73 %求和的上下限為36到36 H=besselh(n-37),2,K*R); %hankel函數(shù) T=(i-(n-37)*exp(i*(n-37)*w(m);% 分子 s(n)=T/H; %對(duì)單個(gè)角度的求和元素endJz(m)=2*sum(s)/(K*Z*pi*R); %單個(gè)角度的電流密度endfigure(1)plot(w*180/pi,abs(Jz),-) %畫(huà)圖，在0360度grid ontitle(TM

20、波入射金屬圓柱體面電流密度分布圖（近似解析解）)xlabel(degrees)ylabel(current dense)Iz=sum(h*Jz) %計(jì)算總電流5.3 TM波入射金屬橢圓柱體的散射function win1(NU) %計(jì)算等效表面電流分布，散射截面分布,NU為波長(zhǎng)%計(jì)算金屬橢圓柱體的散射，橢圓方程：(x/a)2+(y/b)2=1.入射波為T(mén)M波，%電場(chǎng)只有Ez方向入射，用電場(chǎng)積分方程（EFIE）syms x a b; %定義變量Z=377; %特性阻抗K=2*pi/NU; %波數(shù)m=1000; %剖分個(gè)數(shù)P=zeros(m,m); %生成Pmn矩陣框架h=2*pi/m; %剖分角

21、度步長(zhǎng)Q=(0:2*pi/m:2*pi*(m-1)/m); %角度細(xì)分為m等分fai=0; %入射方向和X軸的夾角，可設(shè)定a=NU/4; %橢圓半長(zhǎng)軸b=NU; %橢圓半短軸，當(dāng)a=b時(shí)就是圓柱體散射情形fun=inline(sqrt(a)2*(sin(x).2)+(b)2*(cos(x).2),x,a,b); %橢圓線(xiàn)長(zhǎng)for n=1:m %通過(guò)數(shù)值積分計(jì)算橢圓每個(gè)剖分的線(xiàn)長(zhǎng) x1=(n-1)*h;x2=n*h; %積分上下限delta(n)=quadl(fun,x1,x2,a,b); %積分endX=a*cos(Q+h/2); %對(duì)每個(gè)剖分段，計(jì)算其角度中點(diǎn)的x和y值Y=b*sin(Q+h

22、/2);dx=zeros(m,m); %生成矩陣 dy=zeros(m,m);Del=zeros(m,m);for n=1:m %對(duì)每個(gè)x和y分別計(jì)算(x-xm)2和(y-ym)2dx(n,:)=(X-X(n).2;dy(n,:)=(Y-Y(n).2;Del(n,:)=delta;endI=eye(m);d=dx+dy+I; %(x-xm)2+(y-ym)2,對(duì)角線(xiàn)元素置為1，可以直接用hankel函數(shù)運(yùn)算x=K*sqrt(d); H=besselh(0,2,x); %hankel functionP=Z*K*Del.*H/4; %得到P矩陣Pnn=Z*K*delta.*(1-i*2*log1

23、0(1.78107*K*delta/(4*2.71828)/pi)/4; %計(jì)算對(duì)角線(xiàn)元素 for n=1:m P(n,n)=Pnn(n); %P矩陣對(duì)角元素賦值，endb=exp(-i*K*(X*cos(fai)+Y*sin(fai); %b矩陣,行，要轉(zhuǎn)化為列，M=b.; %b array是復(fù)數(shù)，不能用b取得它的轉(zhuǎn)置，而必須用b.。j=PM; %求解方程，得到電流密度分布arrayas=(j./M)*Z; %abs(j/H),畫(huà)圖的y軸as=as.; figure(1)w2=(0:m/2-1);plot(w2*2/m,abs(as(m/2:-1:1),.) %畫(huà)圖，電流密度在s=0到s=1

24、的分布圖grid ontitle(TM波入射金屬橢圓體的等效表面電流密度分布圖)xlabel(S)ylabel(abs(J/H),current dense/incident magnetic field);Iz=sum(j.*delta) %計(jì)算總電流w=(0:pi/180:pi*179/180);for n=1:180 %計(jì)算散射截面（回波寬度）g=delta.*exp(i*K*(X*cos(w(n)+Y*sin(w(n); %計(jì)算電磁學(xué)課本公式(2.206)的離散方程 G=g.*j.; %j array是復(fù)數(shù)，不能用j取得它的轉(zhuǎn)置，而必須用j.。T=abs(sum(G);out(n)=K

25、*Z2*T2/4; %對(duì)每個(gè)散射角度計(jì)算其散射截面endW=(0:1:179);figure(2)plot(W,sqrt(out/NU),.) %畫(huà)圖，散射截面，0180度grid ontitle(TM波入射金屬圓柱體的散射截面分布圖)xlabel(Degrees)ylabel(square root of echo width/wavelength)5.4 TE波入射金屬橢圓柱體的散射function TEsyms x L S;Z=377;lamda=3;K=2*pi/lamda;a=0.25*lamda; %a,b為橢圓的長(zhǎng)軸和短軸b=1*lamda;N=720; %剖分?jǐn)?shù)目fai=0;

26、%入射角Q=(0:2*pi/N:2*pi*(N-1)/N);h=2*pi/N;X=a*cos(Q+h/2);Y=b*sin(Q+h/2);h=2*pi/N;P=zeros(N,N);fun=inline(sqrt(L)2*(sin(x).2)+(S)2*(cos(x).2),x,L,S); %橢圓線(xiàn)長(zhǎng)for n=1:N %通過(guò)數(shù)值積分計(jì)算橢圓每個(gè)剖分的線(xiàn)長(zhǎng) x1=(n-1)*h;x2=n*h; %積分上下限d=quadl(fun,x1,x2,a,b);%積分delta(n)=d;endfor n=1:N sinn=sign(cos(Q(n)+h/2)*sqrt(1/(a*tan(Q(n)+h/2)/b)2+1); cosn=-1*sign(sin(Q(n)+h/2)*sqrt(1/(1+(b2/(a*tan(Q(n)+h/2)2); Xn=a*cos(Q(n)+h/2);Yn=b*sin(Q(n)+h/2); for m=1:N Xm=a*cos(Q(m)+h/2);Ym=b*sin(Q(m)+h/