基于矩量法的二维金属体散射(内含matlab程序)

基于矩量法的二维金属体散射计算1问题的描述本题是用矩量法计算二维金属圆柱体的散射场，如图所示为一圆柱体和一个椭圆柱的截面，为了计算简单，选入射波为垂直z轴入射的TM或TE平面波yizEizE2x22矩量法求解过程2.1电场积分方程2.1.1问题的分析由麦克斯韦方程组HjE(1)JEjH(2)可得电场积分方程为''20')()(4)(dsKHJKZExz(3)表示在圆柱表面的面电流在远处产生的总场。设入射场为Eiz，散射场为Esz，由金属表面的边界条件szizzEEE=0(4)得''20')()(4)(dlKHJKZECziz(5)2.1.2离散化设入射波为)sincos(yxjkizeE,将散射体截面C分为N份△Cn，用点匹配法对上述积分式子进行离散化，即基函数可取上在其它Cnf10)((6)可得下列离散方程：[P]{J}={b}(7)xy其中：dtyyxxKHKZPmmCmnn))()((42220(8))sincos(imimyxikmeb(9)当m≠n时，))((42220mnmnnmnyyxxKHCKZP(10)当m=n时解析积分为eCKjCZKPnnnn4lg214(11)其中=1.781,e=2.7182.1.3方程组的求解可用LU分解求解方程组，即P=LU，其中P为可逆矩阵，L为上三角矩阵，U为下三角矩阵，则可利用这两个基本的三角矩阵进行求解J，求出J之后，就可求散射场dleyxJKZFFyxjkzCs))sin()cos((),()((12))4/3(81)(KjeKF(13)与二维场中的散射截面2))sin()cos((2),(4)(dleyxJKZyxjKCz(14)2.1.4输出结果的验证此散射问题也可用模式展开法进行求解，可用此结果对本问题进行验证。所得J为nnjnnzKRHejREJ202(15)2.2磁场积分方程对于TE波垂直与z方向入射时的金属体的散射。对于一般的TE波而言只有EEHyxz和,,场分量，电流密度方程只有横向分量。则MFIE为：yAxAJJHxyStSAzttttincz)()()((16)其中')2(0'')(41)()()(dtkRjttttAHJt(17)2'2')()()()(tytytxtxRyy'dltnxx)(sin)(cos)(tytxtt(18)其中t表示边界上的一点，)(t是'dl和X的夹角。根根据前面的过程，圆柱边界分成N分。等效电流密度可以近似为一些脉冲函数的迭加：NnnntttpjJ1)()((19)其中其它上在如果01)(nnCttp(20)则得到HjZincznnmn(21)矩阵非对角元')2(1'''')())()(cos)()((sin4dtkRtyttxtjkmmmCmmmnHRyRxZn(22)2'2')()(tyytxxRmm(23)在nC上认为是常量故)()cos(sin)2(1mmmnmnmnmnkRynxHRyRxZ(24)对角元21Zmn(25)其中22nmnmmnyyxxR回波宽度的近似公式为：21)sincos()sin(4)(NnyxjknnnTEnneCjk(26)3计算机数值实验及分析本论文通过数值计算验证前面理论分析的结果，并对数值计算结果进行分析。分别以金属圆柱体和金属椭圆体为计算例子，做数值实验和分析。所使用的计算机程序是商业软件MATLAB6.5，数值实验在本人机子（celeron41.8GCPU128M内存）,操作系统是windowsxp。3.1二维金属圆柱体的散射基于上面的分析，考虑垂直z方向入射的横向磁波（TM）,离散方程为（7）,编程的基本思路是对(10)式和（11）式编程实现，得出[p]矩阵，再由(9)式得出{b}列，用MATLAB6.5软件上的线性方程组直接求解法求解出{J}。散射截面（回波宽度）可以通过（14）式离散计算出来。计算例子是一个z方向均匀且无限长的金属圆柱，半径为1.5米（5.0R），金属圆柱中心和z轴重合，入射波为z方向极化，幅值为1，从负x轴方向垂直z轴入射的TM平面波，工作频率为100MHz，波长3米。由于是金属体且z方向均匀，可以只考虑对垂直z轴截面的圆周进行剖分并计算。下图给出了720个剖分下电流密度分布的计算结果与近似解析解的比较，其中近似解析解是根据《导波理论》书上3.48式(本文的（15）式)在n=-36到n=36下计算出来的。计算时入射角取为0度。其中x轴为)/(，y轴为电流密度，由图可见，电流密度分布和近似解析解无论幅度相位之间都有着非常好的吻合。计算所得的总等效电流Iz＝-0.0079+0.0083i,而在剖分精度为180时，计算所得的总等效电流Iz=-0.0084+0.0083i。而解析解的总电流Iz=-0.0077+0.0083i，可见随着剖分精度的增加，计算结果收敛于解析解。图1（a）EFIE35.0R，剖分精度720图1（b）EFIE35.0R近似解析解下图给出的是回波宽度的分布图1（c）EFIE35.0R剖分精度720其中x轴为)/(，y轴为dB//2，据个人粗略分析应该基本符合事实。由于没能得到回波宽度的解析解，没能作进一步的分析比较。下面给出入射角为90度，半径2.0R，而其他条件不变的情况下，所得的计算结果。（d）EFIE390,2.0iR剖分精度720由此可见，相对前面那种情况，入射角变化90度，等效电流密度分布也相应有90度的相移，回波宽度的幅度减小了很多，但大体的形状保持不变。这时候的总电流Iz=0.0067+0.0028i。3.2TM波入射金属椭圆柱的散射对于二维金属椭圆柱体的散射这种情况，由于圆柱体是椭圆柱的特殊情况，所以解题的基本思路基本一样，就是对每个剖分步长用数值积分得到，这样有利于得到精确的计算结果。实践的过程也证明了这一点，当每个剖分步长用两个剖分点的直线距离来近似的话，带来很大的误差，而用数值积分得到的结果和解析解很好地吻合。计算例子是一个z方向均匀且无限长的椭圆柱，长轴22b，短轴2/2a，即金属圆柱中心和z轴重合，即椭圆方程为12222byax。入射波为z方向极化，幅值为1，从负x轴方向垂直z轴入射的TM平面波，即入射角为0度。工作频率为100MHz，波长3米。由于是金属体且z方向均匀，可以只考虑对垂直z轴截面的椭圆周进行剖分并计算。图3（a）(b)给出了1000个剖分和2000个剖分下的电流密度分布的计算结果，图3(c)给出解析解以作为比较，而图4(a)(b)给出了上述剖分精度下的回波宽度的计算结果，作为比较图4(c)给出回波宽度的解析解。其中解析解来自参考书《计算电磁场的矩量法》。为了方便与参考书中的解析解比较，x和y轴的参量都相应做了变化。下图中的x轴S为角度的归一化，左端为S=0,右端为S=1。Y轴是iHJ图3(a)EFIE322b2/2a剖分精度1000图3(b)EFIE322b2/2a剖分精度2000图3(c)EFIE322b2/2a剖分精度2700图3(d)EFIE安得列解由上图可见计算结果和参考书提供的解析解能够有很好的吻合，而且随着剖分加细，结果更趋接近于解析解。由于本人机子配置较低，难以对更多剖分点的情况进行运算，但可以预见随着剖分点的增加，计算结果与解析解更好地吻合。但随着剖分数N的增大，计算方法所用的近似不能收敛于解析解，这是因为nm时的，mnP在极限时是不正确的。证明了计算方法的正确性。也可以看到要是用磁场积分方程（MFIE）可以得到更好的解。下图是回波宽度（散射截面）的方向图，其中x轴是角度，y轴是/。图4（a）EFIE322b2/2a剖分精度1000图4(b)EFIE322b2/2a剖分精度2000图4（c）EFIE回波宽度安得列解比较上图可得，计算结果和解析解几乎完全一致，可以注意到即使电流密度分布显著不同，当这两种情况得到的结果却是几乎完全一致的。这是因为/是电流J的连续性泛函，因此J在精确值附近的大小变化是不敏感的。3.3TE波入射金属椭圆柱的散射与上面同样条件，把入射波换为TE波，理论分析可见前面的2.2部分。基于上面的分析，考虑垂直z方向入射的横向电波（TE）,离散方程为（21）,编程的基本思路是对(24)式编程实现，得出[Z]矩阵，再由(17)式得出{H}列，用MATLAB6.5软件上的线性方程组直接求解法求解出{J}。散射截面（回波宽度）可以通过（26）式离散计算出来。图5给出剖分为720份时的计算结果，并给出相应的解析解，以资比较。解析解来自参考书《计算电磁场的矩量法》。图5(a)MFIE322b2/2a剖分精度720图5（b）MFIE安得列解可见，计算结果与解析解大体上符合，但还是存在较大的差别，原因估计是剖分精度不够，还有数值计算[P]矩阵时引进了近似，由于时间仓促，没能在离散化方程时考虑更好的近似方法，这有待于进一步的探讨和研究。图6是TE波入射金属椭圆柱的回波宽度（散射截面）的方向图，其中x轴是角度，y轴是/。可见计算结果和解析解很好地符合，可见虽然电流分布计算结果和所给的解析解有较大的误差，但回波宽度的计算结果和解析解确几乎完全一致，这是因为/是电流J的连续性泛函，因此J在精确值附近的大小变化是不敏感的。图6(a)MFIE322b2/2a剖分精度720图6（b）MFIE安得列解4存在问题和心得存在的问题之一：没考虑到内谐振问题，进一步的工作应该把电场积分方程（EFIE）和磁场积分方程（MFIE）组成联合积分方程（CFIE），来解决这个问题。存在的问题之二：使用的软件MATLAB6.5虽然功能强大，但是运算效率不高，需要占用的内存大和运行时间较长，而比不上用C语言或FORTRAN语言编写的程序效率高。没有利用到课本所介绍的快速多极子技术。存在的问题之三：没有实现在TM波入射情形下用磁场积分方程（MFIE）计算散射场，而据理论分析，用MFIE应该能够得到更好的条件数，计算的结果也能更好地收敛于解析解。这有待于工作的进一步深入。存在问题之四：限于作者知识和经验的不足，没能有更为深厚的理论认识做指导，对结果的分析未免有失偏颇。进一步的工作应该朝着三维散射，介质体散射的方向进行。心得：本文凝结着本组成员的心血，期间经历几多挫折，幸好一一克服了，要说有什么心得的话，第一要对电磁理论有深刻的认识和理解，第二要有深厚的数学基础，对电磁积分方程的性能有深入的理解。第三要熟悉熟练掌握编程语言MATLAB，这是实现的工具。鸣谢：第一应该感谢盛新庆老师的悉心指导，第二感激本组成员的通力合作和不懈的努力。5参考文献[1]盛新庆.计算电磁学要论.北京：科学出版社，2004[2]RogerF.Harrington.FieldComputionbyMomentMethods.NewYork:TheMacmillanCompany,1968[3]AndrewF.Peterson,ScottL.RayandRajMittra.ComputationalMethodsforElectromagnetics.NewYork:IEEEPRESS,1998[4]张志涌等.精MATLAB6.5版.北京：北京航空航天大学出版社，20036程序附录由于作者已经对程序做了很好的注释，应该具备matlab基础的人一般都能够看懂，所以不准备做再多的解释说明。只简单附在后面，以供参考：5.1金属圆柱体散射的程序:functi