10.8线馈矩形微带天线的分析*、**10.8.1三维有限差分法对线馈矩形微带天线的分析**摘要：本文使用三维FDTD 算法实现文献《Application of the three_Dimensional Method to the analysis if Planar Microtrip Circuits 》IEEE trans. On MTT 1990 38(7)的一个矩形微带贴片天线的S11参数的计算。

采用MA TLAB 编程完成数值计算，并与文中的结果进行了比较。

（1） 概述文献《Application of the Three_Dimensional Finite Difference Time Domain Method to the Analysis of Planar Microtrip Circuits 》给出了详细的理论分析。

本文主要是从该文出发，采用MA TLAB 程序完成数值计算过程，画出了时间步为200，400，600，800时介质内的电场分布图形。

天线的尺寸如图10.65所示：图10.70 线馈矩形微带天线结构 （2） 理论基础 支配方程：E t H⨯-∇=∂∂μH tE⨯∇=∂∂ε由此推导出有限差分方程： )()(,1,,,,,1,,,,,,2/1,,,2/1,,,n k j i z n k j i z nk j i y n k j i y n k j i x n k j i x E E yt E E z t H H ---+-∆∆--∆∆+=μμ；*由毕战红, 代子为, 韩春元, 白波, 赵洪涛, 路鹏同学完成)()(1,,,,,,,,1,,,,2/1,,,2/1,,,nk j i x n k j i x n k j i z n k j i z n k j i y n k j i y E E z t E E x t H H ---+-∆∆--∆∆+=μμ； )()(,,1,,,,,1,,,,,2/1,,,2/1,,,nk j i y n k j i y n k j i x n k j i x n k j i z n k j i z E E xt E E y t H H ---+-∆∆--∆∆+=μμ； )()(2/1,,,2/11,,,2/1,,,2/1,1,,,,,1,,,+++++++-∆∆--∆∆+=n k j i y n k j i y n k j i z n k j i z nk j i x n k j i x H H zt H H y t E E εε; )()(2/1,,,2/1,,1,2/1,,,2/1,,1,,,,1,,,+++++++-∆∆--∆∆+=n k j i z n k j i z n k j i x n k j i x nk j i y n k j i y H H xt H H z t E E εε; )()(2/1,,,2/1,1,,2/1,,,2/1,,1,,,,1,,,+++++++-∆∆--∆∆+=n k j i x n k j i x n k j i y n k j i y nk j i z n k j i z H H yt H H x t E E εε （3） 数值计算分析A. 网格划分与时间步确定由于感兴趣的频段范围是DC ——20GHz ，不妨将25GHz 取为频段的上限。

则波长λ的最小值应该是mm f c 12max min ==λ 考虑到∆的取值应该小于等于20min λ所以仅仅从频带的角度考虑，应该有：mm 6.020min max =≤∆λ (10.8.1)z 方向上的介质厚度为0.794mm ，可以将其分为3个网格，z ∆近似有mm z 265.0=∆。

符合max ∆≤∆z 的要求。

Y 方向上的长度为16mm,可以将其分为40个网格或80个网格。

如果分为80个网格，则mm y 2.0=∆，由于Y 方向上的场分布不是我们特别感兴趣的所以不必要将其分的太细，取40个网格就可以。

这样就有mm y 4.0=∆，符合max ∆≤∆y 。

比较困难的是确定x ∆的值，由于在x 方向上有三个尺寸，12.45mm,2.09mm,2.46mm 如果想将每一个尺寸都恰好分为整数个网格数，比较困难。

考虑到天线的尺寸12.45mm 要尽量准确，因此先从这个入手。

文献中给出的389.0=∆x mm ，天线区域分为32个网格，这样有448.1232389.0=⨯，与实际的尺寸有0.002mm 的误差，而微带馈线的宽度为mm 334.26389.0=⨯，误差为-0.126mm,微带馈线的位置为945.15389.0=⨯mm,误差为：-0.115mm.。

可以考虑的另外一种方法：取mm x 2075.0=∆，这样天线区域刚好分为60个网格，没有误差，微带馈线的宽度为12个网格，即mm x 49.212=∆⨯，误差为0.03mm,微带馈线的位置为10个网格，误差为：-0.015mm 。

这样的网格划分可以得到更加精确的模拟，缺点是增加了计算量。

本文采用第一种分法：即将矩形切片尺寸为32x ∆40y ∆⨯3z ∆，总的尺寸为60x ∆⨯100y ∆⨯16z ∆。

确定了z y x ∆∆∆,,以后，可以用稳定性准则确定t ∆。

稳定性条件：2/1222max)111(1-∆+∆+∆≤∆z y x v t 这里t ∆取值为：0.441ps 。

B 源的处理在导带口加强迫激励源，采用了高斯脉冲：220/)(T t t z e E --=.其中0t 为延迟时间，T 高斯脉冲半宽度时间，馈源边界处理为磁壁，T 的取值可由下面公式得出：T=1/f f 为高斯有效频谱的最高频率。

取T=15ps 。

延迟时间0t 取为3T 。

同时为了消除不希望的影响，如虚假反射，该源在存在一定时间后用吸收边界代替（文中取为大于220时间步）。

D 导体的处理在本文中导体看作为无厚度的理导体，在其上的电场切向分量为0。

C 吸收边界采用Mur 一阶吸收边界条件：)(01111n n n n E E yt v y t v E E -∆+∆∆-∆+=++其中 E 0表示网格壁上的切向电场分量，E 1表示为网格内一点的切向电场分量。

D S 参数电场求出后，计算入射波电压iV ][与反射波电压rV ][。

由微波网络理论有i r V S V ]][[][=。

通过傅立叶变化可以求出S 参数：)}({)}({)(t V fft t V fft S i j ij =ω10.8.2三维有限差分法对线馈矩形微带天线的程序与结果*MATLAB 程序：v=3e8;dt=0.441e-12;dx=0.389e-3;dy=0.400e-3;dz=0.265e-3;r=2.2;m=(1+r)/2;A=v*dt/dz; B=v*dt/dy; C=v*dt/dx;D=v*dt/(r*dy);E=v*dt/(r*dz);F=v*dt/(r*dx);G=v*dt/(m*dy);H=v*dt/(m*dz);I=v*dt/(m*dx);J=v*dt/dy;K=v*dt/dz;L=v*dt/dx;a=(v*dt/sqrt(r)-dx)/(v*dt/sqrt(r)+dx);b=(v*dt/sqrt(m)-dx)/(v*dt/sqrt(m)+dx);c=(v*dt-dx)/(v*dt+dx);d=(v*dt/sqrt(r)-dy)/(v*dt/sqrt(r)+dy);e=(v*dt/sqrt(m)-dy)/(v*dt/sqrt(m)+dy);f=(v*dt-dy)/(v*dt+dy);g=(v*dt-dz)/(v*dt+dz);%输入初始值i=2:62;j=2:102;k=2:18;Ex1(i,j,k)=0;Ey1(i,j,k)=0;Ez1(i,j,k)=0;Ex2(i,j,k)=0;Ey2(i,j,k)=0;Ez2(i,j,k)=0;Hx1(i,j,k)=0;Hy1(i,j,k)=0;Hz1(i,j,k)=0;Hx2(i,j,k)=0;Hy2(i,j,k)=0;Hz2(i,j,k)=0;%时间迭代(200,400,600,800)for n=0:200%计算磁场i=2:62;j=2:101;k=2:17;Hx2(i,j,k)=Hx1(i,j,k)+A*(Ey1(i,j,k+1)-Ey1(i,j,k))-B*(Ez1(i,j+1,k) -Ez1(i,j,k));i=2:61;j=2:102;k=2:17;Hy2(i,j,k)=Hy1(i,j,k)+C*(Ez1(i+1,j,k)-Ez1(i,j,k))-A*(Ex1(i,j,k+1) -Ex1(i,j,k));i=2:61;j=2:101;k=2:18;Hz2(i,j,k)=Hz1(i,j,k)+B*(Ex1(i,j+1,k)-Ex1(i,j,k))-C*(Ey1(i+1,j,k) -Ey1(i,j,k));%计算电场i=2:61;j=3:101;k=3:4;%介质层Ex2(i,j,k)=Ex1(i,j,k)+D*(Hz2(i,j,k)-Hz2(i,j-1,k))-E*(Hy2(i,j,k)-H y2(i,j,k-1));i=3:61;j=2:101;k=3:4;Ey2(i,j,k)=Ey1(i,j,k)+E*(Hx2(i,j,k)-Hx2(i,j,k-1))-F*(Hz2(i,j,k)-Hz2(i-1,j,k));i=3:61;j=3:101;k=2:4;Ez2(i,j,k)=Ez1(i,j,k)+F*(Hy2(i,j,k)-Hy2(i-1,j,k))-D*(Hx2(i,j,k)-H x2(i,j-1,k));i=2:61;j=3:101;k=5;%交界面Ex2(i,j,k)=Ex1(i,j,k)+G*(Hz2(i,j,k)-Hz2(i,j-1,k))-H*(Hy2(i,j,k)-H y2(i,j,k-1));i=3:61;j=2:101;k=5;Ey2(i,j,k)=Ey1(i,j,k)+H*(Hx2(i,j,k)-Hx2(i,j,k-1))-I*(Hz2(i,j,k)-H z2(i-1,j,k));i=2:61;j=3:101;k=6:17;%空气层Ex2(i,j,k)=Ex1(i,j,k)+J*(Hz2(i,j,k)-Hz2(i,j-1,k))-K*(Hy2(i,j,k)-H y2(i,j,k-1));i=3:61;j=2:101;k=6:17;Ey2(i,j,k)=Ey1(i,j,k)+K*(Hx2(i,j,k)-Hx2(i,j,k-1))-L*(Hz2(i,j,k)-H z2(i-1,j,k));i=3:61;j=3:101;k=5:17;Ez2(i,j,k)=Ez1(i,j,k)+L*(Hy2(i,j,k)-Hy2(i-1,j,k))-J*(Hx2(i,j,k)-H x2(i,j-1,k));%边界1--接地板i=2:61;j=2:102;Ex2(i,j,2)=0;i=2:62;j=2:101;Ey2(i,j,2)=0;%边界2—微带线贴片i=21:27;j=2:52;Ex2(i,j,5)=0;Ey2(i,j,5)=0;%----矩形切片i=16:48;j=52:92;Ex2(i,j,5)=0;Ey2(i,j,5)=0;%边界3--左侧吸收边界j=2:101;k=3:4;Ey2(2,j,k)=Ey1(3,j,k)+a*(Ey2(3,j,k)-Ey1(2,j,k));j=3:101;k=2:4;Ez2(2,j,k)=Ez1(3,j,k)+a*(Ez2(3,j,k)-Ez1(2,j,k));j=2:101;Ey2(2,j,5)=Ey1(3,j,5)+b*(Ey2(3,j,5)-Ey1(2,j,5));j=2:101;k=6:17;Ey2(2,j,k)=Ey1(3,j,k)+c*(Ey2(3,j,k)-Ey1(2,j,k));j=3:101;k=5:17;Ez2(2,j,k)=Ez1(3,j,k)+c*(Ez2(3,j,k)-Ez1(2,j,k));%边界4--右侧吸收边界j=2:101;k=3:4;Ey2(62,j,k)=Ey1(61,j,k)+a*(Ey2(61,j,k)-Ey1(62,j,k));j=3:101;k=2:4;Ez2(62,j,k)=Ez1(61,j,k)+a*(Ez2(61,j,k)-Ez1(62,j,k));j=2:101;Ey2(62,j,5)=Ey1(61,j,5)+b*(Ey2(61,j,5)-Ey1(62,j,5));j=2:101;k=6:17;Ey2(62,j,k)=Ey1(61,j,k)+c*(Ey2(61,j,k)-Ey1(62,j,k));j=3:101;k=5:17;Ez2(62,j,k)=Ez1(61,j,k)+c*(Ez2(61,j,k)-Ez1(62,j,k));%边界5--后侧吸收边界i=2:61;k=3:4;Ex2(i,102,k)=Ex1(i,101,k)+d*(Ex2(i,101,k)-Ex1(i,102,k)); i=3:61;k=2:4;Ez2(i,102,k)=Ez1(i,101,k)+d*(Ez2(i,101,k)-Ez1(i,102,k)); i=2:61;Ex2(i,102,5)=Ex1(i,101,5)+e*(Ex2(i,101,5)-Ex1(i,102,5)); i=2:61;k=6:17;Ex2(i,102,k)=Ex1(i,101,k)+f*(Ex2(i,101,k)-Ex1(i,102,k)); i=3:61;k=5:17;Ez2(i,102,k)=Ez1(i,101,k)+f*(Ez2(i,101,k)-Ez1(i,102,k));%边界6--上侧吸收边界i=2:61;j=3:101;Ex2(i,j,18)=Ex1(i,j,17)+g*(Ex2(i,j,17)-Ex1(i,j,18));i=3:61;j=2:101;Ey2(i,j,18)=Ey1(i,j,17)+g*(Ey2(i,j,17)-Ey1(i,j,18));%边界7--z方向的棱k=2:17;Ez2(2,2,k)=Ez2(2,3,k)+Ez2(3,2,k)-Ez2(3,3,k);Ez2(62,2,k)=Ez2(62,3,k)+Ez2(61,2,k)-Ez2(61,3,k);Ez2(2,102,k)=Ez2(2,101,k)+Ez2(3,102,k)-Ez2(3,101,k);Ez2(62,102,k)=Ez2(62,101,k)+Ez2(61,102,k)-Ez2(61,101,k);%边界8--x方向的棱i=2:61;Ex2(i,2,18)=Ex2(i,3,18)+Ex2(i,2,17)-Ex2(i,3,17);Ex2(i,102,18)=Ex2(i,101,18)+Ex2(i,102,17)-Ex2(i,101,17);%边界9--y方向的棱j=2:101;Ey2(2,j,18)=Ey2(3,j,18)+Ey2(2,j,17)-Ey2(3,j,17);Ey2(62,j,18)=Ey2(61,j,18)+Ey2(62,j,17)-Ey2(61,j,17);%边界9--源平面if (n>220),i=2:61;k=3:4;Ex2(i,2,k)=Ex1(i,3,k)+d*(Ex2(i,3,k)-Ex1(i,2,k));i=3:61;k=2:4;Ez2(i,2,k)=Ez1(i,3,k)+d*(Ez2(i,3,k)-Ez1(i,2,k));i=2:61;Ex2(i,2,5)=Ex1(i,3,5)+e*(Ex2(i,3,5)-Ex1(i,2,5));i=2:61;k=6:17;Ex2(i,2,k)=Ex1(i,3,k)+f*(Ex2(i,3,k)-Ex1(i,2,k));i=3:61;k=5:17;Ez2(i,2,k)=Ez1(i,3,k)+f*(Ez2(i,3,k)-Ez1(i,2,k));elsei=2:61;k=3:4;Ex2(i,2,k)=Ex1(i,2,k)+D*(2*Hz2(i,2,k))-E*(Hy2(i,2,k)-Hy2(i,2,k-1) );i=2:61;Ex2(i,2,5)=Ex1(i,2,5)+G*(2*Hz2(i,2,5))-H*(Hy2(i,2,5)-Hy2(i,2,4)); i=2:61;k=6:17;Ex2(i,2,k)=Ex1(i,2,k)+J*(2*Hz2(i,2,k))-K*(Hy2(i,2,k)-Hy2(i,2,k-1) );i=3:61;k=2:4;Ez2(i,2,k)=Ez1(i,2,k)+F*(Hy2(i,2,k)-Hy2(i-1,2,k))-D*(2*Hx2(i,2,k) );i=3:61;k=5:17;Ez2(i,2,k)=Ez1(i,2,k)+L*(Hy2(i,2,k)-Hy2(i-1,2,k))-J*(2*Hx2(i,2,k) );i=21:27;k=2:4;Ez2(i,2,k)=exp(-(0.441*n-45)^2/225);endEin(n+1)=Ez2(23,40,5);Esc(n+1)=Ez1(23,40,5)-Ez2(23,40,5);%交换不同时刻的场值Hx1=Hx2;Hy1=Hy2;Hz1=Hz2;Ex1=Ex2;Ey1=Ey2;Ez1=Ez2;endfigure;i=2:62;j=2:102;mesh(Ez1(i,j,3));%%%%%%%%计算s11参数(计算8000步)N=8000;M=0:N;A=fft(Ein);B=fft(Esc);df=1/N/dt;f=M*df;A=abs(A);B=abs(B);s11=20*(log10(A(1:71))-log10(B(1:71))); plot(f(1:71),s11);结果如图10.71到图10.75所图10.71 200时间步图10.72400时间步图10.73600时间步图10.74800时间步图10.75 散射参数S11与文中给出的场分布图形很吻合，其S11参数图形也基本一致。