对称屏蔽带状线——有限差分法

《对称屏蔽带状线——有限差分法》由会员分享，可在线阅读，更多相关《对称屏蔽带状线——有限差分法（17页珍藏版）》请在装配图网上搜索。

1、对称屏蔽带状线的TM模研究,二维有限差分法对称屏蔽带状线的TM模的理论分析1.1 该电磁场边值问题的赫姆霍兹方程及边界条件本文研究的屏蔽带状线结构如图1所示。接地外导体为方形，边长为2b = 5 mm，导带位于中心，平行于上下壁，到左右壁的距离为。令p = a / b,则导带宽度I = 2b(l-p), p的取值范围在0.3到0.7之间，所以导带宽度相对于外导体边长的变化率/ / 2b = 1-p的取值范围也在0.3到0.7之间。图1 对称屏蔽带状线对于TM模，纵向磁场Hz = 0,用纵向电场Ez来描述这个电磁场边值问题，相应的赫姆霍兹方程为(1.1)其中， kc 为对应模式的截止波数。在外导

2、体和导带表面，切向电场分量为零，边界条件为二 0(1.2)z x= bz y = b.二 0y=0,l 2xl 21.2 色散方程及横向场分量的计算公式色散方程为(1.3)求出kc后就可画出两个最低阶TM模的色散特性曲线。求出Ez后可根据以下公式求横向电磁场分量：E = jP. k2 E dxx t czE = jP /k2 -dE dyy. c zH = jgj k2 -dE dyxcz 1H =-js/k2 -dE dxy cz其中，“是对应模式的相位常数。二维有限差分原理lirj- D23hmhzohi1h耳(ij)4边界f图 2 矩形网格的格点配置设y在一个边界C限定的二维场域D内满足

3、亥姆霍兹方程(2.1)6 26 2)+ k2 p 0 .6x2 6y2c 丿首先需将场域D离散化。通常采用比较有规律的分布方式，这样在每个离散点上可得出相同形式的 差分方程，有效地提高解题速度。现以矩形网格的格点配置，导出亥姆霍兹方程的差分形式。设场域D内部某格点0附近的格点分布如图2所示，取步长hi均不相等的情况进行讨论。以附、P、陶、03、P4分别代表在格点0、1、2、3、4处的申的函数值。根据中心差分原理，格点 0 处的二阶偏导的差分形式为V 2申02丿06 2申2丿 0h (申-申)+ h (申-申)沁 2- i01 30h h (h + h )1313h (申-申)+ h (申-申)

4、 沁 2亠20240h h (h + h )2424(2.2)当h = h = h = h = h时，方程(2.1)可差分化为1234(2.3)c0屮+屮 +屮+屮 -=_(hk )2屮12340一般地，可用格点的角标将(2.3)式写为p +p +p+p-4pi+1,ji, j+1i-1,ji,j-1i, j= -(hk)2 屮ci , j(2.4)在波导壁及边界C处，若网格线恰与边界重合，则对任一边界格点b而言，当分析TM模时，相应 的离散化差分格式为申二，b丘C。三、 数值建模3.1 网格划分如图3所示，将外导体包围的场域划分成NxN个正方形网格，网格步长为h二2bN。为使导带正好位于网格

5、点上， N 应取偶数。导带宽度方向的取样点数为(2M +1)，其中M = (b-a)/h ,表示取整。导带中点的坐标为(nj2+1，M 2+1)，最左侧点的坐标为(n2 +1 - M, M 2+1),最右侧点的坐标为(M 2+1 + M,N 2 +】)。3.2差分方程组的导出假设网格划分得足够细，并引入边界条件。现以第2行格点为例进行说明：设九=-(hk)2,2, j3, j,j 二 2,3，,N(3.1)N,j对第 2 行的每个格点都可以写出如(2.4)式的方程一共有N -1个。将方程组写成矩阵的形式,就可以得到(3.2)其中，D是N-1阶方阵，其表示为-41-4-41-4而 I 是 N -

6、1 阶单位方阵。按照同样的方法，依次对第3 行,第4行，第N/2-1行上的格点列出一系列的矩阵方程:I 申 +Dp +I 申=Xp2 3 43I 申 +Dp +Ip =Xp3 454(3.4)I 叫 22+DPN21+I人2=21但是，当运算到第N/2行格点时，由于第N2+1行的格点上存在导带结构(引入0边界条件)，对应的方程组会发生变化，不再与之前的形式相同，具体变为19+ D9 + P9=尢屮N 21N 2N 2+1N 2(3.5)其中，尸是N-1阶对角方阵，P = diag(1.1 00 11),对角线上0元素的个数为 2M+1 。当运算到第 N/2+1 行格点时，由于该行格点包含有导带

7、结构，其场值不再满足方程(2.4)，而是直接为 0。因此，矩阵方程会发生较大的变化，具体形式变为刊 N 2 + 凶 N - + 刊 N -二“ NN 2+1N 2+1(3.6)其中，Q是N-1阶方阵，其表示为41 44004(3.7)414相当于在矩阵D的基础上,将第(i, i+1)和第(i + 1,i)个元素(N/2 M 1 i N/2 + M )由 1 替换成 0，一共发生了 2(2M + 2)次替换，但q矩阵仍然满足对称性。而对于远离导带区域的格点,方程的形式与第2, 3行格点满足的方程形式类似，这里不再重复叙述。将所有的列向量力出，,匕拼成一个长的列向量，再将得到的N个矩阵方程写成方程

8、组的形式，得到K二九(3.8)其中，K为系数矩阵DIIDIK=I D P(3.9)PQPPDIIDIID矩阵K的特征值就是九=-(hk匕,截止波数k对应截止频率f ,即一个TM特征模式，如此c c c就能找出两个最低阶的TM模。特征值九对应的特征向量就是该模式的纵向电场E，根据(1.4)式 z便可以求得横截面上的电/磁场分布。在本问题中，矩阵D, I P, Q均是稀疏对称矩阵，所以矩阵K 是高阶稀疏对称矩阵。在MATLAB中，可利用稀疏矩阵的赋值方式对矩阵K进行存储，进而求解 出最小的两个特征值及其特征向量。这样做能大大降低内存的使用率，同时还能提高计算速度。 3.3双重迭代法求解对于低阶的T

9、M模式，还可采用迭代法1进行求解。求解亥姆霍兹方程与泊松方程不同，因为场 值和特征值-(hk )2均未知，需采用双重迭代法进行求解。首先，除了对各网格点赋予初值，也 c给 k 一个非零初值。然后，按下式对场值进行超松弛迭代：c申n+1 =申ni, ji, j+w14-(kh)2cn + 申 n + 申 n+1 + 申 n+1 )申 ni+1,j i, j+1 i-1,j i, j-1 i, j(3.10)其中，W是迭代加速收敛因子。本问题是正方形场域中的第一类齐次边值问题，最佳的W值可选为w021 + sin (兀 N )(3.11)在经过几次对申 的迭代后，可按下式对特征值-(hk)2进行迭

10、代i , j c工申 Cp +申 +申 +申知 )ASi, j i, ji, j+1i1,ji, j1i, ji, j(hk)2 = -irjcp 2 ASi,j i,ji,j其中，AS代表第(i, j)个格点所占的面积，对于内点AS = h2。接着，再利用求得的新k值代入i, ji,jc(3.10)式中，迭代求得新的,如此重复这种双重迭代过程，直到场值收敛到应有的精度为止，程i,j序框图如图 4 所示。经过多次尝试后发现，若对待求的网格点作简单、常规的初始值设定，利用双 重迭代法求解时，最终都将收敛至最低阶的TM模式。因此，双重迭代法只适用于求解最低阶的TM 模，很难系统地求解更高阶的模式。

11、若要求得更高阶的TM模式，还是需要对系数矩阵K进行特征 值和特征向量的求解。因此，本文将采用两种计算方法(双重迭代法和特征值法)计算最低阶的TM 模式，采用特征值法计算次低阶的TM模式。开始否MadtN4)=0?是是辙出结果按超松弛迭代法进行 次这代，柚厂设呂场埋内炉的初11和的韧值蛤定边界衆件对廐的值进行迭代迭诧次数计數Z结束图 4 双重迭代法程序框图四、 结果分析41导带宽度对两个最低阶TM模的影响网格数为30x30时，两个最低阶TM模式的截止频率随p的变化关系如表1所示，与文献2中 的计算结果基本相等。随着p的增大，広减小，/2不变；即随着导带宽度的减小，最低阶TM模的 截止频率降低，次

12、低阶TM模的截止频率不变。表1中列举了用两种方法计算出的広，可以看出二 者几乎相等，这也从侧面验证了二维有限差分法计算程序的正确性。表1截止频率随p的变化（N = 30）计算方法p0.30.40.50.60.7双重迭代法fc133.421533.188932.920331.893231.1024特征值法fc133.420833.188232.919531.892231.1014fc233.465833.465833.465833.465833.46584.2两个最低阶TM模的色散特性和横截面场分布由表1可知，p取不同值时，c2不变，Cl的变化也不大，由式（1.3）可知色散特性曲线几乎不变。只需

13、画出p取某个特定值时两个TM模的色散特性曲线。图5是p = 0.5时两个最低阶TM模的色散特性曲线。图5 p = 0.5时两个最低阶TM模的色散特性曲线不管p取何值，两个最低阶TM模的横向场分布不变，只是幅度略有不同。这里给出p = 0.5、工作频率为40 GHz时对称屏蔽带状线横截面的场分布，如图6和7所示，网格数为30x30。仔细观察这两个模式的电场和磁场分布，可以发现，最低阶 TM 模的电场和磁场的相位在导带所在的平面 发生 180的跃变，导带宽度决定了相位跃变区域的大小。而次低阶 TM 模的场分布与方波导中的 TM12模相同，不受导带影响，因此，当导带宽度减小即p增大时，最低阶TM模的

14、截止频率逐渐下 降，接近方波导中TM11的截止频率（约为21.2 GHz），而次低阶TM模截止频率不变。-SCS-50图6最低阶TM模的横截面场分布图:（a）电场;（b）磁场.图 7 次低阶 TM 模的横截面场分布图: (a) 电场 ; (b) 磁场 . 4.3 网格粗细对结果的影响如表2所示，取p = 0.5,研究网格粗细对结果的影响。表2网格粗细对结果的影响(p = 0.5)计算方法N2030405060双重迭代法4132.748832.920332.710832.822232.6842特征值法fc132.747832.919532.710232.821732.6838fc233.4008

15、33.465833.488533.499133.5048取N = 20至60之间的偶数进行数值实验，并画出fc1和厶随N值变化的曲线图。如图8所示， 随着网格数目的增加，fCi的值呈现出振荡的趋势，但总体趋势是在逐渐减小；而fc2则是缓慢增大, 并趋近于一个稳定值(约为33.51 GHz)。一般来说，当网格足够细时，计算结果应该趋于某个稳定值，fc2的变化趋势很好地符合了这个规律，但的变化趋势不太符合该规律。当取p = 0.5时，令N = 2n(N为偶数)。本程序在进行模拟时，采用的是“四舍五入”的取整方 式，由 1.1 和 3.1 可知，此时 M 的取值为(4.1)n / 2, n e ev

16、en(n +1)/2, n e odd以N的取值为20, 22, 24的情况来进行说明。此时对应的M的取值为5, 6, 6。实际计算时，导 带宽度相对于外导体边长的比值为2M/ N,对应的取值为0.5, 0.55, 0.5。由此可见，当N = 22时， 程序所模拟的导带宽度大于预设宽度，由4.1可知，広随导带宽度的增大而升高；当N = 24时，导 带的相对宽度又再次减小为0.5,厶便随之降低，这也就解释了图8(a)中产生振荡现象的原因。而fc2 几乎不随导带宽度的变化而变化，因此不会产生数值振荡的现象。広和/2的变化趋势也可能与数值 色散有关，即随着数值波模在网格中的传播方向以及离散化情况不同

17、，传播速度将随频率改变。(a)(b)图8 截止频率随网格数目的变化. (af)c1; (b)fc2五、 总结本文以基于二维有限差分法的程序代码为工具研究介绍了对称屏蔽带状线的导带宽度对两个最 低阶TM模截止频率的影响。计算结果表明，随着导带宽度的减小，最低阶TM模的截止频率逐渐 降低，次低阶TM模的截止频率几乎不变。本文还对网格粗细对计算结果的影响做了数值实验，发 现有限差分法的网格尺寸要在计算精度和数值色散之间折中，一般每个波长长度划分10-20个网格 点，此处的波长是指主要频谱分量对应的波长。参考文献1 王秉中，计算电磁学M.北京：科学出版社，2002.2 刘静，赵德双，王秉中.基于二维有

18、限差分法的屏蔽带状线的TM模研究C/全国微波毫米波 会议. 2007.附录%*特征值方法求两个最低阶的 TM 模* clear;clc;%*设置基本常数* c=2.99792458e8;muz=4.0*pi*1.0e-7;epsz=1.0/(c 人2*muz);%speed of light in free space %permeability of free space %permittivity of free space%*设置尺寸和网格参数* b=5e-3;N=30;h=2*b/N;rho=0.5;a=b*rho;M=round(b-a)/h);%外导体边长的一半 %网格步数 %网格

19、步长 %a/b 的取值%导带宽度方向的取样点数为2M+1;%*设置分块系数矩阵* K=cell(N-1);D=zeros(N-1);I=eye(N-1,N-1);P=eye(N-1,N-1);%*建立矩阵 D*for i=1:N-1D(i,i)=-4;endfor i=1:N-2D(i,i+1)=1;D(i+1,i)=1;end%*建立矩阵 P* for i=N/2-M:N/2+MP(i,i)=0;end%*建立矩阵 Q* Q=D;for i=N/2-M-1:N/2+MQ(i,i+1)=0;Q(i+1,i)=0;end %*转换成稀疏矩阵存储* D=sparse(D);I=sparse(I);

20、P=sparse(P);Q=sparse(Q);%*建立 K 矩阵for i=1:N-1for j=1:N-1Ki,j=sparse(N-1,N-1);endendfor i=1:N-1Ki,i=D;endfor i=1:N-2Ki,i+1=I;Ki+1,i=I;endKN/2,N/2=Q;KN/2-1,N/2=P;KN/2+1,N/2=P;KN/2,N/2-1=P;KN/2,N/2+1=P;%*求特征值与特征向量 V,D=eigs(K,2,sm);lamda=diag(D);kc=sqrt(-lamda)/h; fc=c*kc/(2*pi*1e9);Ez2=reshape(V(:,2),N-

21、1,N-1); u2=zeros(N+1);f=40e9;omega=2*pi*f;k=2*pi*f/c;beta2=sqrt(k 人2-kc(2)人2);K=cell2mat(K);%将原胞数组转换成普通数组%求最小的两个特征值及对应的特征向量%将K的特征向量转换为格点的场值%工作频率%工作波数%相位常数，kc(2)为次低阶TM模的截止波数 %若要计算最低阶TM模式的横向场分布，只需将相应的被注释语句取消注释即可% Ez1=reshape(V(:,1),N-1,N-1);% beta1=sqrt(k 人2-kc(l)人2);% u1=zeros(N+1);%*计算横截面电/磁场分布Ex=ze

22、ros(N+l);Ey=zeros(N+1);Hx=zeros(N+1);Hy=zeros(N+1);for i=2:Nfor j=2:Nu2(i,j)=Ez2(i-1,j-1);%引入边界格点上的0 场值%u1(i,j)=Ez1(i-1,j-1);endend for i=2:Nfor j=2:NEx(i,j)=-lj*beta2/kc(2)人2*(u2(i+l,j)-u2(i-l,j)/(2*h);Ey(i,j)=-lj*beta2/kc(2)人2*(u2(i,j+l)-u2(i,j-l)/(2*h);Hx(i,j)=lj*omega*epsz/kc(2)人2*(u2(i,j+l)-u2(

23、i,j-l)/(2*h);Hy(i,j)=-lj*omega*epsz/kc(2)人2*(u2(i+l,j)-u2(i-l,j)/(2*h);%Ex(i,j)=-lj*betal/kc(l)人2*(ul(i+l,j)-ul(i-l,j)/(2*h);%Ey(i,j)=-lj*betal/kc(l)人2*(ul(i,j+l)-ul(i,j-l)/(2*h);%Hx(i,j)=lj*omega*epsz/kc(l)人2*(ul(i,j+l)-ul(i,j-l)/(2*h);%Hy(i,j)=-lj*omega*epsz/kc(l)人2*(ul(i+l,j)-ul(i-l,j)/(2*h);end

24、end %*计算边界点处的横向场值Ex(l,:)=0;Ex(N+l,:)=0;Ey(:,l)=0;Ey(:,N+l)=0;%切向电场为零Hy(l,:)=0;Hy(N+l,:)=0;Hx(:,l)=0;Hx(:,N+l)=0;%法向磁场为零%*用平均值的方法计算边界处的法向电场和切向电场*Ey(l,2:N)=2*Ey(2,2:N)-Ey(3,2:N);Ey(N+l,2:N)=2*Ey(N,2:N)-Ey(N-l,2:N);Ex(2:N,l)=2*Ex(2:N,2)-Ex(2:N,3);Ex(2:N,N+l)=2*Ex(2:N,N)-Ex(2:N,N-l);Hx(l,2:N)=2*Hx(2,2:N

25、)-Hx(3,2:N);Hx(N+l,2:N)=2*Hx(N,2:N)-Hx(N-l,2:N);Hy(2:N,l)=2*Hy(2:N,2)-Hy(2:N,3);Hy(2:N,N+l)=2*Hy(2:N,N)-Hy(2:N,N-l);%*最/次低阶 TM 模式横向场分布图* x=le3*linspace(-b,b,N+l);y=le3*linspace(-b,b,N+l);figure(l)quiver(x,y,imag(Ey),imag(Ex),k);xlabel(x/mm),ylabel(y/mm);axis equal;axis(le3*-b b -b b);title( (a)电场);f

26、igure(2) quiver(x,y,imag(Hy),imag(Hx),k); xlabel(x/mm),ylabel(y/mm);axis equal;axis(1e3*-b b -b b);%这里画出的是次低阶的TM模场分布图x=30:0.1:42;yl=sqrt(2*pi*x/c*le9).人2-kc(l)人2);y2=sqrt(2*pi*x/c*le9).人2-kc(2)人2); figure(3)plot(x,real(yl),-.k,x,real(y2),r);xlabel(频率 f(GHz),ylabel(相位常数beta(m人-1); title(两个最低阶TM模色散曲线) legend(最低阶TM模,次低阶TM模);grid on;