电磁场矩量法解金属圆柱导体散射

矩量法解无限长金属圆柱导体散射场摘要矩量法（MOM）在研究电磁散射的问题中十分有效，能帮我们解决很多电磁场的问题。本文主要对电磁场的矩量法进行了学习和研究，并利用Matlab进行编程仿真，对特定问题进行求解，求解出无限长圆柱导体的雷达散射截面，并对结果进行分析。0引言电磁场在我们的生活中扮演者越来越重要的作用。我们手机的信号要靠电磁场来进行传播，很多探测类的工作都要靠电磁场的探测来进行实现。特别是在军事中，探测飞机的雷达散射信号来获取敌方军情和设计己方飞机不被探测到，即检测和设计隐形飞机有着十分重要的意义。所以我们要对物体对电磁场的散射情况进行分析，主要是通过麦克斯韦方程这一基本公式来进行各种计算。而由于实际问题十分复杂，导致解析方法难以实现，所以需要数值计算的方法来进行实现。利用计算机强大的计算能力，可以迅速的算出电磁场对物体的散射情况，再对其进行分析，以得到我们所需要的场，从而探测到物体的形状位置等重要信息。其中重要的一个环节，就是对问题的分析和方法的应用，而矩量法在电磁散射中有着广泛的应用，本文就将对应用矩量法来求解无限长圆柱导体的雷达散射问题。1矩量法概述矩量法（TheMethodofMoments,简称MOM），在天线、微波技术和电磁波散射等方面应用十分广泛。矩量法最早被Richmand和Harrington用于求解电磁场的问题，而后在Harrington的著作中得到了系统的论述，从此成为求解电磁场问题数值解的主要方法，并成功应用于天线和电磁散射问题，至今已有50年的发展历史。矩量法是一种将连续方程离散化为代数方程组的方法。其基本原理是：先选定基函数对未知函数进行近似展开，带入算子方程，再选取适当的权函数，使在加权平均的意义下方程的余量等于零，由此将连续的算子方程转换为代数方程。剩下的问题就是利用计算机进行大量的数字计算。原则上，矩量法可用于求解微分方程和积分方程，但用于微分方程时所得到的的系数矩阵往往是病态的，故在电磁场中主要用于求解积分方程。矩量法是一种严格的数值方法，求解精度高，加之格林函数直接满足辐射条件，无需设置吸收边界条件，因而可以灵活解决边界比较复杂的一些问题，在电磁辐射和散射、天线电流分布、天线设计、微波网络、生物电磁学、辐射效应研究、微带线分析、电磁兼容等方面得到广泛应用。下面本文就将对TM波入射无限长金属圆柱导体的雷达散射截面采用矩量法做进一步的计算分析。2问题分析矩量法把泛函方程转化为矩阵方程，然后通过矩阵方程求解。一般表达式会是如下的泛函方程：（1）φ是一个函数，求解一个未知函数比较困难，我们熟悉的是求解方程中的未知数而不是未知函数，所以我们可以将其展开为级数的形式，如下：（2）这样（1）变为：（3）两边乘权函数并积分得：（4）可以写为：（5）这样可以表示为一个矩阵形式：（6）其中，Amn=wm，L(vn)m,n=1,2,…,N（7）Bm=wm,fm,n=1,2,…,N（8）取基函数为：（9）权函数为：（10）可得：（11）bm=φ(12)这样就可以求出{c}：c=A-1b（13）则各点的场就能求解。已知一个无限长的导体圆柱，有TM波Ei射向圆柱体，求散射后的场分布，只考虑二维情况。图一示意图TM极化波电场满足亥姆霍兹方程：（14）在导体圆柱表面，由边界条件有：（15）（16）还有方程：（17）（18）由上述方程我们可以得到：（19）（20）即可得散射场在圆柱体表面和远处的分布。对其进行数值化：（21）（22）（23）可得：（24）用计算机进行编程运算，即可求得解。3Matlab编程仿真主要流程如下：具体做法可见附录，有详细的文字说明。这里简要说明一下圆柱圆周的分割问题。示意图如图所示：主程序几何分割添矩阵解矩阵方程计算并绘图End图二圆柱划分示意图先将其分为N等分，求出各点的坐标，再取每段的中点，求出新的坐标，避免端点的奇异性：（25）每段长度间隔：（26）参数设置：波长为1，入射波表达式为：（27）其中入射角为π。4仿真结果及分析运行程序得到了导体圆柱表面电流密度分布图（图三），散射场图（图四）和总场图（图五），还有散射动态图，无法展示可运行查看结果。其参数设置都为入射角为π，半径r=λ，图片所显示的分割精度为100次，可自行调节分割次数，次数越多精度越高，但速度也会越慢。图三TM波入射的金属圆柱表面等效电流分布图四TM波入射金属圆柱散射场分布图五TM波入射金属圆柱总场分布图图六课本上的例图对比仿真图与课本上的从图中我们可以看到仿真结果基本正确，与书上的结果大致吻合。电流密度在入射的方向最大，背对入射方向最小，符合基本原理。但是在背对入射场的方向，即角度为0度附近，仿真的等效电流不均匀，与解析解有一定的偏差，可能是精度不够，或者是奇异点问题处理的不够好，有待进一步改进。此题采用Mie级数进行求解，得到解析解为：nnjnnzKRHejREJ202（28）利用matlab对其进行绘图，程序见附录二，得到的等效电流分布如图七所示：图七Mie级数近似解析解与之前采用矩量法得到的等效电流分布图出了刚刚所说的端点处，基本完全吻合，再一次证明了所求解的场分布的正确性。但是当波长变短，即频率增加时，电流分布会变得极其不稳定，但改变剖分精度，又能基本稳定，如图八和图九所示。参数分别表在图下方。图八λ=1/3.5,N=100图九λ=1/3.5,N=500个人认为有可能是剖分精度不同，导致有些点在计算中边为奇异点，导致解发生震荡。其中原理还没研究太明白。5总结通过此次矩量法的研究和学习，成功地仿真出无限长金属圆柱导体对TM波散射情况的分布图，收获很大。在刚开始接触矩量法时，完全没有思路。后来通过大量的文献查阅和学习，逐渐理解了矩量法的基本思路和编程方法。在仿真的过程中，也遇到过各种问题，挑战很大，后来通过反复看文献理解矩量法的思想，学习Matlab的编程方法才慢慢有了一定的体会，最终才得到了较为正确的结果，也有了对矩量法更为深刻的认识。这次的仿真，不仅仅是一次完成作业，更是一次挑战，不仅学到了课本上的知识，还学到了Matlab的编程技巧，收获远不止一次完成任务这么简单。虽然仿真结果基本吻合，但是此次仿真还有很多不足之处。首先是程序计算速度还不够快，未能结合数值方法对其进行加速处理，程序有待进一步优化。其次就是未能对导体内部共振对散射场的影响做出讨论和分析，结果可能不是特别准确。最后，由于时间比较紧，理论推导中的公式比较复杂，所以没有花时间进行公示的编辑，只是在书上进行截图，有些地方可能不清晰，有待改善和提高。最后感谢老师这一学期认真负责的为我们讲解电磁场数值计算的知识，经常在一些我们困惑的地方知道我们，学到了很多，感谢李老师！参考文献[1]JinJ.TheoryandComputationofElectromagneticFields[M]//Theoryandcomputationofelectromagneticfields.Wiley:,2010.[2]佚名.电磁场数值计算[M].高等教育出版社,1996.[3]何红雨.电磁场数值计算法与MATLAB实现[J].2004.[4]JinJ.电磁场矩量法[M].西安电子科技大学出版社,1998.[5]王长清.现代计算电磁学基础[J].2005.[6]陈涌频,孟敏,方宙奇.电磁场数值方法[M].科学出版社,2016.[7]胡光华.计算电磁学要论[J].2008(10):3-4.附录一矩量法程序functionMOMwavelength=1;%设置入射波长rho=1;%导体圆柱半径Z0=377;%真空波阻抗k0=2*pi/wavelength;%波矢b=(Z0*k0/4);%计算要用到的常量N=100;%导体圆柱所分段数Z_mn=zeros(N,N);%定义阻抗矩阵V_m=zeros(N,1);%定义入射场矩阵width=10*rho;%作图区域宽度Nx=100;%作图划分精度Ny=100;%****************************划分圆柱导体*****************************forn=1:N+1%分成N段，求每点的坐标x(n)=rho*cos(2.0*pi*(n-1.0)/N);y(n)=rho*sin(2.0*pi*(n-1.0)/N);endforn=1:N%求两点中点坐标，避免端点奇异性xc(n)=(x(n+1)+x(n))/2;yc(n)=(y(n+1)+y(n))/2;L(n)=sqrt((x(n+1)-x(n))^2+(y(n+1)-y(n))^2);%分段长度end%****************************求解矩阵**************************************form=1:NV_m(m)=Einc(xc(m),yc(m),k0);%不同点的入射场赋值forn=1:N%带入分析后的式子if(m==n)%考虑奇异点情况Z_mn(m,n)=b*L(n)*(1-1j*(2/pi)*(log(0.455*k0*L(n))-1));elser=sqrt((xc(m)-xc(n))^2+(yc(m)-yc(n))^2);Z_mn(m,n)=b*L(n)*besselh(0,2,k0*r);endendendJ=Z_mn\V_m;%求解表面电流分布figure(1)%画出电流密度在0－360度的分布图w=(1:1:m);plot(w*360/m,abs(J))gridon;xlabel('degrees');ylabel('Js/H0');%添加网格线,设置坐标轴title('TM波入射的金属圆柱表面等效电流密度分布图')%***************************计算并画散射场场图******************************A=zeros(Nx,Ny);%定义散射场矩阵B=zeros(Nx,Ny);%定义总场矩阵forix=1:Nx%画一定区域内的场图xo=(2*((ix-1)/(Nx-1))-1)*width;foriy=1:Nyyo=(2*((iy-1)/(Ny-1))-1)*width;forn=1:NR=sqrt((xo-xc(n))^2+(yo-yc(n))^2);A(ix,iy)=A(ix,iy)-b*J(n)*L(n)*besselh(0,2,k0*R);endendendAoutput=abs(A);save('scattered.dat','Aoutput','-ASCII')%保存数据figure(2)%画散射场分布图pcolor(Aoutput');title('TM波入射金属圆柱散射场图');colorbar;shadinginterp;axisequal;xlim([1100]);%***************************计算并画总场图**********************************forix=1:Nxxo=(2*((ix-1)/(Nx-1))-1)*width;foriy=1:Nyyo=(2*((iy-1)/(Ny-1))-1)*width;forn=1:NR=sqrt((xo-xc(n))^2+(yo-yc(n))^2);B(ix,iy)=B(ix,iy)-b*J(n)*L(n)*besselh(0,2,k0*R);endB(ix,iy)=B(ix,iy)+Einc(xo,yo,k0);endendBoutput=abs(B);save('total.dat','Boutput','-ASCII')%保存数据figure(3)%画总场的分布图pcolor(Boutput');title('TM波入射金属圆柱总场图');colorbar;shadinginterp;axisequal;xlim([1100]);%****************************总场动态图*******************************