設計任務采用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計算出滿足條件的是