基于时域有限差分法(FDTD)的矩形谐振腔分析

采用FDTD数值计算的方法来分析理想谐振腔中的场，谐振腔尺寸为25*12.5*60mm填充空气，采用直角坐标系下的场分量迭代公式，激励源采用高斯脉冲源，源的参数根据谐振腔的尺寸来确定。分析时间和空间离散度以及采样点数对分析结果的影响。

二、 方案设计

（1）学习FDTD理论，并推导直角坐标系下maxwell方程的差分方程；

（2）理论学习并推导理想矩形谐振腔中的时谐场，并分析其谐振频率分布； （3）激励源采用高斯脉冲源，导体采用PEC边界，利用FDTD编程求解谐振腔内的场分量；

（4）对谐振腔内部分点处的采样数据进行频谱分析，提取其谐振频率分布，并与理论对比，并分析时间和空间离散度以及采样点数对分析结果的影响。

三、 设计原理

3.1时域有限差分法

FDTD(finite diference time domain)方法属于全波分析法， 它是Yee在1966年所提出的数值方法“ ，其原理是将麦克斯韦方程式中两个微分形式的旋度方程式以中心差分式做离散化。求解过程由递推完成，尤其适合计算机编程实现。 3.1.1有限差分法

有限差分法是用变量离散的、含有有限个未知数的差分方程近似的代替连续变量的微分方程，即构造合理的差分格式，使其解能保持原问题的主要性质，并有相当高的精确度。假设f(x),为x的连续函数，在x轴上每隔h距离取一点，其中任意某一点用xi表示，则叫做f(x)在xi点的中心差分。在时域有限差分法中正是用中心差商代替微商，同时用Max-well方程组建立差分方程。

3.1.2 Yee’s差分算法

E,H 场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理

HDEtt(i+1,j,k+1)Ez HxEyHzExExEyHyEzEyEz(i+1,j+1,k+1)ExEBH tt --（1）

如图3-1所示，Yee单元有以下特

点：

（1）E与H分量在空间交叉放

置，相互垂直；

(i,j+1,k+1)xzy(i,j,k)(i,j+1,k)（2）每一坐标平面上的E分量四周由H分量环绕，H分量的四周由E分量环绕；

（3）每一场分量自身相距一个空 间步长，E和H相距半个空间步长 （4）电场取n时刻的值，磁场取n+0.5时刻的值；

（5）电场n+1时刻的值由 n 时刻的值得到，磁场n+0.5时刻的值由n-0.5时刻的值得到；电场n+1时刻的旋度对应(n+1)+0.5时刻的磁场值，磁场n+0.5时刻的旋度对应 (n+0.5)+0.5时刻的电场值； （6）3个空间方向上的时间步长相等，以保证均匀介质中场量的空间变量与时间变量完全对称。

时域Maxwell方程租由两个旋度方程和两个散度方程构成，两个旋度方程是安培环路定律和法拉第电磁感应定律的微分形式，概括了宏观电磁场的基本规律。时域Maxwel方程租在直角坐标中，可以写成6个标量方程

--(2)

--(3)

E、H、、、、磁场强度、介电常数、磁 导率、电导率、导磁率。

以上6个偏微分方程是FDTD算法的基础。1966年，K． S．Yee 提出了一种E、 分量的节点在空间上的交替排列和在 时间上的交替抽样方式，从而可以在时间轴上逐次推进的求解空问电磁场的值。正是靠这种合理的Yee氏网格体系，才成功的创立了FDTD算法。

在这个矩阵差分网格单元中，定义网格单元顶点的空间坐标(x,y,z)=(ix,iy,iz)简写为(i,j,k)其中x，y,z别表示在(x,y,z)坐标方向的网格步长，(i,j,k)为整数。在时间上取时间步长为t,则电场分量在n时刻的时间 =nt；而磁场分量应在与电场分量相差半个时间步长处取样，即磁场的取样点为 te-t= (n-)t。根据时间和空间网格划分的规律,任意～个空间和时间的电场或磁场分量函数可表示为

)=F(ix, jy, kz, nt)=F’’(i,j,k) -- (4)

采用中心差分来替代对时间、空间坐标的微分，得到

-- (5)

可以看出，在任意时间步长上空间任意网格点上的电场值取决于3个因素：①该点在上一时间步长的电场值；②与该电场正交平面上邻近点处在上一时间步长上的磁场值；③介质的电参数和 。磁场值同理。通过这些基本算式，在逐个时间步长对模拟区域各网格点的电、磁场交替进行计算，在执行到适当的时间步数后，即可获得需要的时域数值结果。这种差分格式通常称之为蛙跳格式。

在FDTD算法的显示蛙跳格式中，每一步计算都无需做矩阵求逆运算，避免了矩阵求逆运算带来的许多问题，这是该方法的一个突出优点。

3.2谐振腔谐振频率的理论分析

在谐振腔内部，电磁波频率为驻波分布，即在空间中有固定的波节波腹，

只有一些特定的频率能够建立起稳定的驻波，从而实现谐振，这些频率成为谐振频率。谐振腔中，电磁波频率只能取不连续的离散值。谐振腔内波的波数为

无论是在谐振腔中还是在波导中，本征方程, 都必须成立，否则波动方程没有非零解。在波导中和是由边界条件确定的，因此，对于任意工作频率总存在对应的传播常数，总存在相应的本征解(也就是说电磁场总能在波导中存在)。虽然这个解可能是传输状态也可能处于截止状态。

谐振频率是谐振腔的主要参数之一，矩形谐振腔中TE、TM 模式的谐振频率具有相同的数学表达式：

-- (6)

矩形谐振腔的谐振频率是腔体几何尺寸(宽a,高b,长l）的函数。由于金属腔体的几何尺寸会随环境温度变化而改变，而且，微波、毫米波的中心工作频率很高，所以环境温度的变化将会导致谐振频率绝对值的较大变化

如果激励源的频率高于该模式的截止频率时，该模式就成为传输模式；如果激励源的频率低于该模式的截止频率时，该模式就成为截止模式。对于微波谐振腔而言，由于它比传输线多了纵向边界条件，要想在微波谐振腔中激励起某个模式，激励源不但要与该模式满足奇偶禁戒规则，而且激励源的频率必须等于该模式的谐振频率。

① 波导中，不管微波源的频率与波导的截止频率关系如何，波动方程 总是有解的，这个解可能是截止状态也可能是传输状态。 ②在谐振腔中，如果微波源的频率不等于谐振腔的谐振频率，波动方程就没有非零解。产生这一区别的根本原因是纵向边界条件。

③谐振腔有无穷多个谐振频率，每个谐振频率对应一个或多个模式。如果不同的模式(m，n，p 中任一个或多个参数不同，则电磁场的空间分布就不同。)具有相同的谐振频率，这些模式就称为简并模式。

简并模式的特点就是它们具有相同的谐振频率和不同的电磁场空间分布结构。在一般情况下，我们不希望微波谐振腔工作在简并模式。但在某些特殊场合，可以利用简并模式的特点达到特殊的目的。如TM模，m,n=0,1,2,…;但不能同时为0；TE模，p不能为0

由于TM模和TE模都能存在于矩形波导内，所以TM模和TE模也同样可以存在于矩形谐振腔中。由于谐振腔内不存在唯一的纵方向（传播方向），因此，TM模和TE模的名称不唯一。假设z轴为参考的“传播方向”，由于在z=0和z=l处存在导体壁，电磁波将在期间来回反射形成驻波，所以在 空腔内不可能有波的传播。

3. 3 三维问题（直角坐标系）

先讨论式第一个公式：

HyExHzE --(7)

yzt11在ttn2时间步，对节点 i2,j,k的离散公式为：

εi,j,k1211Exn1i2,j,kExni2,j,kΔt Hinz1212γi,j,k12En1xi12,j,kEi,j,k2nx12,j,kHzΔy1212n12iny12121,j2,k Hyn12i,j,k12HiΔz121,j,k2

-- (8)

上式中的第二项用平均值来替代Exn12i12,j,k是因为离散方程中电场的时间取

样是整数n，磁场的时间取样是n+1/2，所以只能取n及n+1时电场的平均值。实际也证明这个平均值使FDTD算法具有数值稳定性。

1n1整理后，将Exi2,j,k作为未知数，其余作为迭代计算的已知数 1Exn1i2,j,kCAmExni12,j,kCBm 1111nnn2n1111111Hzi2,j2,kHz2i2,j2,kHy2i2,j,k2Hy2i2,j,kΔyΔz12

1式中：mi2,j,k --（9）

mt2m2 CAmt

mmmt1t22mmm1t1m CBm mmmt1t22m同理，式（1-3）中其它两个公式的离散形式为

11n1nEyi,j2,kCAmEyi,j2,kCBm 1111nnnn211111111222Hi,j,kHi,j,kHi,j,kHi,j,kxxzz22222222zx1式中：mi,j2,k --（10） 11Ezn1i,j,k2CAmEzni,j,k2CBm 1111nnnn211111111222 Hi,j,kHi,j,kHi,j,kHi,j,kyyxx22222222xy1式中：mi,j,k2 --（11）

以上三式是电场的时间推进计算公式。

同样，讨论式（3-4）中第一个公式，设观察点x,y,z为Hx的节点，即在时

11刻tnt，对节点 i,j2的离散公式为： ,k2HxCPmHi,j,k Ei,j1,kEi,j,kEi,j,k1Ei,j,k CQm n1212i,j,k12nx121212nz12nz12ny12ny12yz11式中，mi,j2 ----（12） ,k2mmt2mt2 CPm mmtmmm1t22mmmm1t1mCQm mmmt1t22m同理，式（1-4）中其它两个公式的离散形式为

HyCPmHi,j,k Ei,j,k1Ei,j,kEi1,j,kEi,j,k CQm n1212i,j,k12ny1212nx12nx1212nz12nz12zx11式中，mi2 --（13） ,j,k2Hz11,j,kCPmHzi2,j2,k 1111nnnnEyi1,j2,kEyi,j2,kExi2,j1,kExi2,j,k CQm xy1212n12in1211式中，mi2,j2,k --（14）

以上三式是磁场的时间推进计算公式。 时域推进计算框图（交叉半步逐步推进）

若已知t1t0nt时空间各节点处的电场值（赋初值） 计算t2t1t时空间各节点处的磁场值，式（8）—（11） 2计算t1t2

在编程中，为了使电场和磁场有相同的数量级（为减小误差），可对H或E~~进行“归一化”处理，即：用HZ0H取代H，用EE/Z0取代E，式中Z0是自由空间波阻抗。计算结果再分别除以和乘以Z0即可。

可以看出，这种离散方法电场和磁场在时间顺序上交替抽样，抽样间隔相差半个时间步长，使麦克斯韦方程离散后成为显示差分方程，从而可以在时间上迭代求解，不需矩阵求逆。给定初值后，可以逐步推进，求得以后各个时刻点的空间电磁分布。这是FDTD法的最大特点。

t时空间各节点处的电场值，式（8）—（10） 23.4 FDTD法的稳定性及数值色散问题

在FDTD中，时间增量t和空间增量x、y、z之间不是相互独立的，它们的取值必须满足一定的关系，以避免数值结果的不稳定——表现为随着时间步数的增加，计算结果发散。造成解不稳定的因素有多种：

① 计算机在计算过程中，原始数据可能有误差，如系数阵建立过程中产生

的误差，而每次运算由于只能保留有限位数而又产生误差，误差的积累有可能淹没真正解，使计算结果不可靠，即不稳定； ② 计算方法不合适； ③ 离散间隔不当等。

为了确定数值解稳定的条件，有许多推导方式，结论相同. 1、时间步长稳定性要求（推导略）

一般情况下 tT1 （15） T为波动周期。

2、时间步长与空间步长的关系

三维 t1111xyz222 （16）

在非均匀区域，v取最大值。真空中v=c(光速)。若是立方体Yee元胞，那么txyz，

13 （ 17）。若是正方形Yee元胞（二维），xy，

t1那么t12

（18） 。若是线段等分Yee元胞（一维），那么

 （19）

若电磁波所在空间的介质特性与频率有关，则电磁波的传播速度也将是频率的函数，这种现象称色散。而FDTD方程是原Maxwell方程的一种近似，所以当计算机对电磁波在空间的传播进行模拟时，在非色散空问中也会出现色散现象， 且电磁波的相速度随波长、传播方向及变量离散化的情况发生变化，这种非物理性的色散称为数值色散。数值色散会导致脉冲波形的破坏，出现人为的各向异性及虚假的折射现象。数值色散是由于近似差商替代连续微商引起的，这种影响可 以通过减小离散化过程所取空间和时间步长而无限减小，但计算空间的总网格数目的增加也会相应增加对计算机存储空间和计算时间的要求，所以要根据实际条件来选择合适的步长。而数值色散在FDTD法中是不可避免的。

为了减小数值色散，除了满足式（16）外，可取比式（15）更高的要求

tT （20） 12并令电磁波沿网格的对角线方向传播，这时有kxkykzk，可以得到理想的色散关系。

3.5 建模及原程序

% %矩形腔 宽度为25mm,高度为12.5mm,谐振腔长度60mm，介质为空气 % % 计算得谐振腔的谐振频率为13.42Ghz clear;clc;

c=3e8; muz=4.0*pi*1.0e-7; epsz=1.0/(c^2*muz); %每小格长,三个方向的网格长均设为dx dx=0.00125;

X=25e-3;Y=12.5e-3;Z=60e-3;

ie=X/dx; je=Y/dx; ke=Z/dx; ib=ie+1; jb=je+1; kb=ke+1; %源及观察的位置

is=ie/2+1; js=je/2+1; ks=ke/2+1; %时间步长

dt=dx/(2*c); f0=13.42e9; %总步数

nmax=[800,1000，3000]; , %% 选择f=1; nmax=nmax(f); %脉冲的系数

rtau=0.3/f0;%决定脉宽的时间长度 tau=rtau/dt;%归dt化 ndelay=tau*3;%时延 %介质的参数

%选择腔内任意一点进行验证 i=10;j=5;k=10; %介质的参数

epsr=1.0; sig=0.0; mur=1.0; sigm=0.0; %迭代系数

ca=(1.0-(dt*sig)/(2.0*epsz*epsr))/(1.0+(dt*sig)/(2.0*epsz*epsr)); cb=(dt/epsz/epsr/dx)/(1.0+(dt*sig)/(2.0*epsz*epsr));

da=(1.0-(dt*sigm)/(2.0*muz*mur))/(1.0+(dt*sigm)/(2.0*muz*mur));db=(dt/muz/mur/dx)/(1.0+(dt*sigm)/(2.0*muz*mur)); %初始化

ex=zeros(ie,jb,kb);ey=zeros(ib,je,kb);ez=zeros(ib,jb,ke); hx=zeros(ib,je,ke);hy=zeros(ie,jb,ke);hz=zeros(ie,je,kb); for n=1:nmax

ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));

ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));

ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke)); %高斯线源脉冲,电流源沿着z方向激励

ez(is,js,1:ke)=ez(is,js,1:ke)+8*exp(-((n-ndelay)^2/tau^2));

hx(2:ie,1:je,1:ke)=da*hx(2:ie,1:je,1:ke)+db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));

hy(1:ie,2:je,1:ke)=da*hy(1:ie,2:je,1:ke)+db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));

hz(1:ie,1:je,2:ke)=da*hz(1:ie,1:je,2:ke)+db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke)); %将某点处的值取出来

o(1,n)=ex(i,j,k);o(2,n)=ey(i,j,k);o(3,n)=ez(i,j,k); o(4,n)=hx(i,j,k);o(5,n)=hy(i,j,k);o(6,n)=hz(i,j,k);

end

%选择定点进行谐振频率观察

fs=1/dt; N=floor((nmax-1)/2); pl=(0:(N-1))*fs/nmax ; xs={'ex','ey','ez','hx','hy','hz'}; for I=1:6

I=1:6

F(I,:)=abs(fft(o(I,:),nmax)); subplot(3,2,I)

plot(pl/1e9,F(I,1:N));

set(gca,'xlim',[0,50]) xlabel('频率 Ghz') legend(xs{I}); grid on

end

四、 结果及分析

矩形波导谐振腔是由两端短路的一段金属波导构成，矩形波导参数----宽度a=25 mm，高度b=12.5 mm，长度l= 60 mm。由微波理论可知，当b（1）决定网格单元的尺寸和时间步长取Δx=Δy=Δz=Δ=1.25mm，则空间网格

数为50×25×120，取Δt=Δ/(2c)=2.0833ps。

（2）设置边界条件对于矩形波导谐振腔，腔体的6个面都是金属，并且均为理

想导体，即腔内导体边界上的所有切向电场分量为0，所有法向磁场分量为0。矩形谐振腔的谐振波长为： ，则TE101模的谐振波长为=22.3606mm，主模谐振频率为=13.42Ghz

（3）设置激励源高斯脉冲的表达式为: 为了在谐振腔

中激励起模，并且抑制其他高次模，选择线源脉冲，使之在腔内xy平面中

心处沿y轴方向分布，并选择合适的和τ值。经过反复试验，取=138.462ps，τ=46.154ps。

1ex0-1 0400ez2000 00.1hy0.050 0 1ey0-1 00.2hx0.10 01hz0-1 0 103020频率 Ghz4050 103020频率 Ghz4050 103020频率 Ghz4050 103020频率 Ghz4050 103020频率 Ghz4050103020频率 Ghz4050 图 4-1 分析时长为（dt*nmax=2.0833ps*800）注意对应时，在谐

振腔体中

上图显示了在谐振腔内部由高斯脉冲源激励起由谐振腔自身参数选择出的谐振频率分量，（12.5，6.25，12.5）mm处各个长分量所包含的频率分量与理论值对比，发现在谐振频率13.2Ghz处出产生了最大响应。波形图中的其他波峰是由于高斯激励脉冲元不是理想脉冲所引起的，对本次验证试验没有太大影响。因为不是绝对理想情况，与计算所得的13.42数据误差相差不大，即得到验证。

4） 时间和空间离散度以及采样点数对分析结果的影响

时间离散度不但应满足数值色散的最低要求，而且要与空间补偿满足稳定关系，在这次试验中，根据Yee元胞模型，在三个坐标方向均采用等步长dx=1.25mm，根据3-12,3-15计算出满足条件的是时间步长dt=2.0833ps；总的时间步数为nmax，设计了三个备选值（800,1000,3000）。

采样有效时间越长，即在相同的采样频率（时间步长dt的倒数）下，采样点数nmax增加。那么频域中一个周期内频率的间隔fo（1/dt/nmax）减小，

对两个频率相近序列的频谱峰值分辨能力增强，谱峰更尖锐更细，频率分辨力越好

1ex0-1 02000ez10000 00.2hy0.10 0 1ey0-1 01hx0.50 01hz0-1 0 102030频率 Ghz4050 102030频率 Ghz4050 102030频率 Ghz4050 102030频率 Ghz4050 102030频率 Ghz4050102030频率 Ghz4050 图 4-2 分析时长为（dt*nmax=2.0833ps*3000）时，在谐振腔体中 （12.5，6.25，12．5）mm处各场分量激起的归一化的幅频响应图

五、 结论

麦克斯韦方程组是支配宏观电磁现象的一组基本方程，这组方程有微分和积分形式；FDTD方法就是由微分形式的麦克斯韦旋度方程出发进行差分离散从而得到一组时域推进公式。使用matlab对已知尺寸填充空气理想导体边界的矩形谐振腔进行三维的maxwell旋度方程分析，使用高斯脉冲激励源激励，验证了物体的电尺寸决定了物体的电磁谐振频率（固有频率）和腔内存在的主模，均与外激励无关。另外，在离散空间中，当时间和空间步长确定后，随着采样点数的增加，幅频响应途中频率分辨率变高，这是在数字信号处理中经常应用到的概念。以上结合FDTD和MATLAB对矩形波导谐振腔做了仿真分析，程序简洁明了，运行效率也较高。FDTD法在电磁场数值分析方面有很大的优越性，而matlab处理数据和图形十分好用，二者结合可将算法迅速程序化，并能获得很好的数据处理结果。

六、 参考文献

1.《电磁场与电磁波》， Matlab相关书籍等 2. 葛德彪，闫玉波，《电磁波时域有限差分方法》，西安电子科技大学出版

3. 基于FDTD法的电磁仿真的优化算法 谭文泉 (南京工业大学信息科学与工程

学院，江苏南京210009)

4. 部分程序参考了Susan C. Hagnes 公布在网上的源程序Susan C. Hagness Department of Electrical and Computer EngineeringUniversity of Wisconsin-Madison