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

[1]Jin J. Theory and Computation of Electromagnetic Fields[M]// Theory and computation of electromagnetic fields. Wiley :, 2010. [2] 佚名. 电磁场数值计算[M]. 高等教育出版社, 1996. [3] 何红雨. 电磁场数值计算法与 MATLAB 实现[J]. 2004. [4] Jin J. 电磁场矩量法[M]. 西安电子科技大学出版社, 1998. [5] 王长清. 现代计算电磁学基础[J]. 2005. [6] 陈涌频, 孟敏, 方宙奇. 电磁场数值方法[M]. 科学出版社, 2016. [7] 胡光华. 计算电磁学要论[J]. 2008(10):3-4.
用于天线和电磁散射问题， 至今已有 50 年的发展历史。 矩量法是一种将连续方程离散化为代数方程组的方法。 其基本原理是：先选定基函数对未知函数进行近似展开， 带入算子方程，再选取适当的权函数，使在加权平均的 意义下方程的余量等于零，由此将连续的算子方程转换 为代数方程。剩下的问题就是利用计算机进行大量的数 字计算。原则上，矩量法可用于求解微分方程和积分方 程，但用于微分方程时所得到的的系数矩阵往往是病态 的，故在电磁场中主要用于求解积分方程。矩量法是一 种严格的数值方法，求解精度高，加之格林函数直接满 足辐射条件，无需设置吸收边界条件，因而可以灵活解 决边界比较复杂的一些问题，在电磁辐射和散射、天线 电流分布、天线设计、微波网络、生物电磁学、辐射效 应研究、微带线分析、电磁兼容等方面得到广泛应用。 下面本文就将对 TM 波入射无限长金属圆柱导体的雷达 散射截面采用矩量法做进一步的计算分析。 2 问题分析 矩量法把泛函方程转化为矩阵方程，然后通过矩 阵方程求解。一般表达式会是如下的泛函方程： （1）
附录一 矩量法程序 function MOM wavelength=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; %****************************划分圆柱导体***************************** for n=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); end for n=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 %****************************求解矩阵************************************** for m=1:N V_m(m)=Einc(xc(m),yc(m),k0); %不同点的入射场赋值 for n=1:N %带入分析后的式子 if (m==n) %考虑奇异点情况 Z_mn(m,n)=b*L(n)*(1-1j*(2/pi)*(log(0.455*k0*L(n))-1)); else r=sqrt((xc(m)-xc(n))^2+(yc(m)-yc(n))^2); Z_mn(m,n)=b*L(n)*besselh(0,2,k0*r); end end end J=Z_mn\V_m; %求解表面电流分布 figure(1) % 画出电流密度在 0－360 度的分布图 w=(1:1:m); plot(w*360/m,abs(J)) grid on;xlabel('degrees');ylabel('Js/H0'); %添加网格线,设置坐标轴 title('TM 波入射的金属圆柱表面等效电流密度分布图') %***************************计算并画散射场场图****************************** A=zeros(Nx,Ny); %定义散射场矩阵 B=zeros(Nx,Ny); %定义总场矩阵 for ix=1:Nx %画一定区域内的场图 xo=(2*((ix-1)/(Nx-1))-1)*width; for iy=1:Ny
φ 是一个函数，求解一个未知函数比较困难，我们熟悉
（14） 在导体圆柱表面，由边界条件有： （15）
的是求解方程中的未知数而不是未知函数，所以我们可 以将其展开为级数的形式，如下： （2） 这样 （1） 变为： 两边乘权函数并积分得： （4） 可以写为： （5） 这样可以表示为一个矩阵形式： （6） 其中， Amn=<wm， L(vn)> Bm=<wm,f> 取基函数为： （9） 权函数为： （ 10 ） 可得： （11） bm=φ 这样就可以求出{c}： c=A-1b 则各点的场就能求解。 已知一个无限长的导体圆柱， 有 TM 波 Ei 射向圆柱 体，求散射后的场分布，只考虑二维情况。 （13） (12) m,n=1,2,…,N m,n=1,2,…,N （7） （8） （3）
End 图一 示意图 TM 极化波电场满足亥姆霍兹方程： 具体做法可见附录，有详细的文字说明。这里简要说 明一下圆柱圆周的分割问题。示意图如图所示：
图三 TM 波入射的金属圆柱表面等效电流分布
图二 圆柱划分示意图 先将其分为 N 等分，求出各点的坐标，再取每段的中 点，求出新的坐标，避免端点的奇异性： （25） 每段长度间隔： （26） 参数设置： 波长为 1，入射波表达式为： （27） 其中入射角为π。 4 仿真结果及分析 运行程序得到了导体圆柱表面电流密度分布图 （图三） ，散射场图（图四）和总场图（图五） ，还有 散射动态图，无法展示可运行查看结果。其参数设置 都为入射角为 π，半径 r=λ，图片所显示的分割精度为 100 次，可自行调节分割次数，次数越多精度越高，但 速度也会越慢。 图五 TM 波入射金属圆柱总场分布图 图四 TM 波入射金属圆柱散射场分布
图七 Mie 级数近似解析解 与之前采用矩量法得到的等效电流分布图出了刚刚所 说的端点处，基本完全吻合，再一次证明了所求解的 场分布的正确性。 图六 课本上的例图 对比仿真图与课本上的从图中我们可以看到 仿真结果基本正确，与书上的结果大致吻合。电流密 度在入射的方向最大，背对入射方向最小，符合基本 原理。但是在背对入射场的方向，即角度为 0 度附 近，仿真的等效电流不均匀，与解析解有一定的偏 差，可能是精度不够，或者是奇异点问题处理的不够 好，有待进一步改进。 此题采用 Mie 级数进行求解，得到解析解为： 但是当波长变短，即频率增加时，电流分布会变 得极其不稳定，但改变剖分精度，又能基本稳定，如 图八和图九所示。参数分别表在图下方。
（16） 还有方程： （17） （18） 由上述方程我们可以得到： （19） （20） 即可得散射场在圆柱体表面和远处的分布。 对其进行数值化： （21） （22） （23） 可得：
（24） 用计算机进行编程运算，即可求得解。 3 Matlab 编程仿真 主要流程如下： 主程序
几何分割 添矩阵 解矩阵方程 计算并绘图
Jz
2 E0 j n e jn 2 KR R.5, N=100
利用 matlab 对其进行绘图，程序见附录二，得到的等 效电流分布如图七所示：
图九 λ=1/3.5, N=500
个人认为有可能是剖分精度不同，导致有些点在 计算中边为奇异点，导致解发生震荡。其中原理还没 研究太明白。 5 总结
yo=(2*((iy-1)/(Ny-1))-1)*width; for n=1:N R=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); end end end Aoutput=abs(A); save('scattered.dat','Aoutput','-ASCII') %保存数据 figure(2) % 画散射场分布图 pcolor(Aoutput');title('TM 波入射金属圆柱散射场图'); colorbar;shading interp;axis equal;xlim([1 100]); %***************************计算并画总场图********************************** for ix=1:Nx xo=(2*((ix-1)/(Nx-1))-1)*width; for iy=1:Ny yo=(2*((iy-1)/(Ny-1))-1)*width; for n=1:N R=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); end B(ix,iy) = B(ix,iy)+Einc(xo,yo,k0); end end Boutput=abs(B); save('total.dat','Boutput','-ASCII') %保存数据 figure(3) % 画总场的分布图 pcolor(Boutput');title('TM 波入射金属圆柱总场图'); colorbar;shading interp;axis equal;xlim([1 100]); %****************************总场动态图************************************* figure(4) for k=1:90 % 设置运行时间为 90*0.1=9s phi=2*pi*k/30; C=real(B*exp(1j*phi)); pcolor(C');title('TM 波入射金属圆柱总场动态图'); colorbar;shading interp; axis equal; xlim([1 100]); pause(0.1); %0.1s 刷新一次 end %******************************设置 TM 波入射角度***************************** function E=Einc(x,y,k0) phi=pi; % 入射角度设为π E=exp(1j*k0*(x*cos(phi)+y*sin(phi))); %入射场 TM 波表达式，幅度为 1