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

设计任务采用FDTD数值计算的方法来分析理想谐振腔中的场，谐振腔尺寸为25*12.5*60mm填充空气，采用直角坐标系下的场分量迭代公式，鼓励源采用高斯脉冲源，源的参数根据谐振腔的尺寸来确定。分析时间和空间离散度以及采样点数对分析结果的影响。方案设计〔1〕学习FDTD理论，并推导直角坐标系下maxwell方程的差分方程；〔2〕理论学习并推导理想矩形谐振腔中的时谐场，并分析其谐振频率分布；〔3〕鼓励源采用高斯脉冲源，导体采用PEC边界，利用FDTD编程求解谐振腔内的场分量；〔4〕对谐振腔内局部点处的采样数据进行频谱分析，提取其谐振频率分布，并与理论比照，并分析时间和空间离散度以及采样点数对分析结果的影响。设计原理3.1时域有限差分法FDTD(finitediferencetimedomain)方法属于全波分析法，它是Yee在1966年所提出的数值方法“，其原理是将麦克斯韦方程式中两个微分形式的旋度方程式以中心差分式做离散化。求解过程由递推完成，尤其适合计算机编程实现。3.1.1有限差分法有限差分法是用变量离散的、含有有限个未知数的差分方程近似的代替连续变量的微分方程，即构造合理的差分格式，使其解能保持原问题的主要性质，并有相当高的精确度。假设f(x),为x的连续函数，在x轴上每隔h距离取一点，其中任意某一点用xi表示，则叫做f(x)在xi点的中心差分。在时域有限差分法中正是用中心差商代替微商，同时用Max-well方程组建立差分方程。3.1.2Yee’s差分算法场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理--〔1〕如图3-1所示，Yee单元有以下特点：〔1〕与分量在空间交叉放置，相互垂直；〔2〕每一坐标平面上的分量四周由分量环绕，分量的四周由分量环绕；〔3〕每一场分量自身相距一个空间步长，和相距半个空间步长〔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三维问题〔直角坐标系〕先讨论式第一个公式：--(7)在时间步，对节点的离散公式为：--(8)上式中的第二项用平均值来替代是因为离散方程中电场的时间取样是整数n，磁场的时间取样是n+1/2，所以只能取n及n+1时电场的平均值。实际也证明这个平均值使FDTD算法具有数值稳定性。整理后，将作为未知数，其余作为迭代计算的数式中：--〔9〕同理，式〔1-3〕中其它两个公式的离散形式为式中：--〔10〕式中：--〔11〕以上三式是电场的时间推进计算公式。同样，讨论式〔3-4〕中第一个公式，设观察点为的节点，即在时刻，对节点的离散公式为：式中，----〔12〕同理，式〔1-4〕中其它两个公式的离散形式为式中，--〔13〕式中，--〔14〕以上三式是磁场的时间推进计算公式。时域推进计算框图〔交叉半步逐步推进〕假设时空间各节点处的电场值假设时空间各节点处的电场值〔赋初值〕计算时空间各节点处的磁场值，式〔8〕—〔11〕计算时空间各节点处的电场值，式〔8〕—〔10〕在编程中，为了使电场和磁场有相同的数量级〔为减小误差〕，可对H或E进行“归一化〞处理，即：用取代，用取代，式中是自由空间波阻抗。计算结果再分别除以和乘以即可。可以看出，这种离散方法电场和磁场在时间顺序上交替抽样，抽样间隔相差半个时间步长，使麦克斯韦方程离散后成为显示差分方程，从而可以在时间上迭代求解，不需矩阵求逆。给定初值后，可以逐步推进，求得以后各个时刻点的空间电磁分布。这是FDTD法的最大特点。3.4FDTD法的稳定性及数值色散问题在FDTD中，时间增量和空间增量、、之间不是相互独立的，它们的取值必须满足一定的关系，以防止数值结果的不稳定——表现为随着时间步数的增加，计算结果发散。造成解不稳定的因素有多种：计算机在计算过程中，原始数据可能有误差，如系数阵建立过程中产生的误差，而每次运算由于只能保存有限位数而又产生误差，误差的积累有可能淹没真正解，使计算结果不可靠，即不稳定；计算方法不适宜；离散间隔不当等。为了确定数值解稳定的条件，有许多推导方式，结论相同.1、时间步长稳定性要求〔推导略〕一般情况下〔15〕T为波动周期。2、时间步长与空间步长的关系三维〔16〕在非均匀区域，v取最大值。真空中v=c(光速)。假设是立方体Yee元胞，，那么〔17〕。假设是正方形Yee元胞〔二维〕，，那么〔18〕。假设是线段等分Yee元胞〔一维〕，那么〔19〕假设电磁波所在空间的介质特性与频率有关，则电磁波的传播速度也将是频率的函数，这种现象称色散。而FDTD方程是原Maxwell方程的一种近似，所以当计算机对电磁波在空间的传播进行模拟时，在非色散空问中也会出现色散现象，且电磁波的相速度随波长、传播方向及变量离散化的情况发生变化，这种非物理性的色散称为数值色散。数值色散会导致脉冲波形的破坏，出现人为的各向异性及虚假的折射现象。数值色散是由于近似差商替代连续微商引起的，这种影响可以通过减小离散化过程所取空间和时间步长而无限减小，但计算空间的总网格数目的增加也会相应增加对计算机存储空间和计算时间的要求，所以要根据实际条件来选择适宜的步长。而数值色散在FDTD法中是不可防止的。为了减小数值色散，除了满足式〔16〕外，可取比式〔15〕更高的要求〔20〕并令电磁波沿网格的对角线方向传播，这时有，可以得到理想的色散关系。3.5建模及原程序%%矩形腔宽度为25mm,高度为12.5mm,谐振腔长度60mm，介质为空气%%计算得谐振腔的谐振频率为13.42Ghzclear;clc;c=3e8;muz=4.0*pi*1.0e-7;epsz=1.0/(c^2*muz);%每小格长,三个方向的网格长均设为dxdx=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);forn=1:nmaxex(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'};forI=1:6I=1:6F(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});gridonend结果及分析矩形波导谐振腔是由两端短路的一段金属波导构成，矩形波导参数----宽度a=25mm，高度b=12.5mm，长度l=60mm。由微波理论可知，当b<a<l时，TE101模的谐振波长最长，是矩形波导谐振腔的主模。TE波〔即横电波〕是在波的传播方向上有磁场分量，但没有电场分量，即电场垂直于传播方向。TE101模在矩形腔的3个方向都不传输能量，呈驻波分布，并且由边界条件决定：电场分量只有Ez分量，在腔体中央最强；磁场有Hx和Hy两个分量，在腔壁附近最强，腔体中央为0。由于矩形波导谐振腔有着广泛的应用，并且大多工作在主模状态，这里选择微分高斯脉冲以鼓励起模。〔1〕决定网格单元的尺寸和时间步长取Δx=Δy=Δz=Δ=1.25mm，则空间网格数为50×25×120，取Δt=Δ/(2c)=2.0833ps。〔2〕设置边界条件对于矩形波导谐振腔，腔体的6个面都是金属，并且均为理想导体，即腔内导体边界上的所有切向电场分量为0，所有法向磁场分量为0。矩形谐振腔的谐振波长为：QUOTE，则TE101模的谐振波长为=22.3606mm，主模谐振频率为=13.42Ghz〔3〕设置鼓励源高斯脉冲的表达式为:QUOTEQUOTE为了在谐振腔中鼓励起模，并且抑制其他高次模，选择线源脉冲，使之在腔内xy平面中心处沿y轴方向分布，并选择适宜的和τ值。经过反复试验，取=138.462ps，τ=46.154ps。图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计算出满足条件的是