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

設計任務采用FDTD數(shù)值計算的方法來分析理想諧振腔中的場，諧振腔尺寸為25*12.5*60mm填充空氣，采用直角坐標系下的場分量迭代公式，鼓勵源采用高斯脈沖源，源的參數(shù)根據(jù)諧振腔的尺寸來確定。分析時間和空間離散度以及采樣點數(shù)對分析結果的影響。方案設計〔1〕學習FDTD理論，并推導直角坐標系下maxwell方程的差分方程；〔2〕理論學習并推導理想矩形諧振腔中的時諧場，并分析其諧振頻率分布；〔3〕鼓勵源采用高斯脈沖源，導體采用PEC邊界，利用FDTD編程求解諧振腔內(nèi)的場分量；〔4〕對諧振腔內(nèi)局部點處的采樣數(shù)據(jù)進行頻譜分析，提取其諧振頻率分布，并與理論比照，并分析時間和空間離散度以及采樣點數(shù)對分析結果的影響。設計原理3.1時域有限差分法FDTD(finitediferencetimedomain)方法屬于全波分析法，它是Yee在1966年所提出的數(shù)值方法“，其原理是將麥克斯韋方程式中兩個微分形式的旋度方程式以中心差分式做離散化。求解過程由遞推完成，尤其適合計算機編程實現(xiàn)。3.1.1有限差分法有限差分法是用變量離散的、含有有限個未知數(shù)的差分方程近似的代替連續(xù)變量的微分方程，即構造合理的差分格式，使其解能保持原問題的主要性質(zhì)，并有相當高的精確度。假設f(x),為x的連續(xù)函數(shù)，在x軸上每隔h距離取一點，其中任意某一點用xi表示，則叫做f(x)在xi點的中心差分。在時域有限差分法中正是用中心差商代替微商，同時用Max-well方程組建立差分方程。3.1.2Yee’s差分算法場分量取樣節(jié)點在空間和時間上采取交替排布，利用電生磁，磁生電的原理--〔1〕如圖3-1所示，Yee單元有以下特點：〔1〕與分量在空間交叉放置，相互垂直；〔2〕每一坐標平面上的分量四周由分量環(huán)繞，分量的四周由分量環(huán)繞；〔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個空間方向上的時間步長相等，以保證均勻介質(zhì)中場量的空間變量與時間變量完全對稱。時域Maxwell方程租由兩個旋度方程和兩個散度方程構成，兩個旋度方程是安培環(huán)路定律和法拉第電磁感應定律的微分形式，概括了宏觀電磁場的根本規(guī)律。時域Maxwel方程租在直角坐標中，可以寫成6個標量方程--(2)--(3)E、H、、、、磁場強度、介電常數(shù)、磁導率、電導率、導磁率。以上6個偏微分方程是FDTD算法的根底。1966年，K．S．Yee提出了一種E、分量的節(jié)點在空間上的交替排列和在時間上的交替抽樣方式，從而可以在時間軸上逐次推進的求解空問電磁場的值。正是靠這種合理的Yee氏網(wǎng)格體系，才成功的創(chuàng)立了FDTD算法。在這個矩陣差分網(wǎng)格單元中，定義網(wǎng)格單元頂點的空間坐標(x,y,z)=(ix,iy,iz)簡寫為(i,j,k)其中x，y,z別表示在(x,y,z)坐標方向的網(wǎng)格步長，(i,j,k)為整數(shù)。在時間上取時間步長為t,則電場分量在n時刻的時間=nt；而磁場分量應在與電場分量相差半個時間步長處取樣，即磁場的取樣點為te-t=(n-)t。根據(jù)時間和空間網(wǎng)格劃分的規(guī)律,任意～個空間和時間的電場或磁場分量函數(shù)可表示為)=F(ix,jy,kz,nt)=F’’(i,j,k)--(4)采用中心差分來替代對時間、空間坐標的微分，得到--(5)可以看出，在任意時間步長上空間任意網(wǎng)格點上的電場值取決于3個因素：①該點在上一時間步長的電場值；②與該電場正交平面上鄰近點處在上一時間步長上的磁場值；③介質(zhì)的電參數(shù)和。磁場值同理。通過這些根本算式，在逐個時間步長對模擬區(qū)域各網(wǎng)格點的電、磁場交替進行計算，在執(zhí)行到適當?shù)臅r間步數(shù)后，即可獲得需要的時域數(shù)值結果。這種差分格式通常稱之為蛙跳格式。在FDTD算法的顯示蛙跳格式中，每一步計算都無需做矩陣求逆運算，防止了矩陣求逆運算帶來的許多問題，這是該方法的一個突出優(yōu)點。3.2諧振腔諧振頻率的理論分析在諧振腔內(nèi)部，電磁波頻率為駐波分布，即在空間中有固定的波節(jié)波腹，只有一些特定的頻率能夠建立起穩(wěn)定的駐波，從而實現(xiàn)諧振，這些頻率成為諧振頻率。諧振腔中，電磁波頻率只能取不連續(xù)的離散值。諧振腔內(nèi)波的波數(shù)為無論是在諧振腔中還是在波導中，本征方程,都必須成立，否則波動方程沒有非零解。在波導中和是由邊界條件確定的，因此，對于任意工作頻率總存在對應的傳播常數(shù)，總存在相應的本征解(也就是說電磁場總能在波導中存在)。雖然這個解可能是傳輸狀態(tài)也可能處于截止狀態(tài)。諧振頻率是諧振腔的主要參數(shù)之一，矩形諧振腔中TE、TM模式的諧振頻率具有相同的數(shù)學表達式：--(6)矩形諧振腔的諧振頻率是腔體幾何尺寸(寬a,高b,長l〕的函數(shù)。由于金屬腔體的幾何尺寸會隨環(huán)境溫度變化而改變，而且，微波、毫米波的中心工作頻率很高，所以環(huán)境溫度的變化將會導致諧振頻率絕對值的較大變化如果鼓勵源的頻率高于該模式的截止頻率時，該模式就成為傳輸模式；如果鼓勵源的頻率低于該模式的截止頻率時，該模式就成為截止模式。對于微波諧振腔而言，由于它比傳輸線多了縱向邊界條件，要想在微波諧振腔中鼓勵起某個模式，鼓勵源不但要與該模式滿足奇偶禁戒規(guī)則，而且鼓勵源的頻率必須等于該模式的諧振頻率。波導中，不管微波源的頻率與波導的截止頻率關系如何，波動方程總是有解的，這個解可能是截止狀態(tài)也可能是傳輸狀態(tài)。②在諧振腔中，如果微波源的頻率不等于諧振腔的諧振頻率，波動方程就沒有非零解。產(chǎn)生這一區(qū)別的根本原因是縱向邊界條件。③諧振腔有無窮多個諧振頻率，每個諧振頻率對應一個或多個模式。如果不同的模式(m，n，p中任一個或多個參數(shù)不同，則電磁場的空間分布就不同。)具有相同的諧振頻率，這些模式就稱為簡并模式。簡并模式的特點就是它們具有相同的諧振頻率和不同的電磁場空間分布結構。在一般情況下，我們不希望微波諧振腔工作在簡并模式。但在某些特殊場合，可以利用簡并模式的特點到達特殊的目的。如TM模，m,n=0,1,2,…;但不能同時為0；TE模，p不能為0由于TM模和TE模都能存在于矩形波導內(nèi)，所以TM模和TE模也同樣可以存在于矩形諧振腔中。由于諧振腔內(nèi)不存在唯一的縱方向〔傳播方向〕，因此，TM模和TE模的名稱不唯一。假設z軸為參考的“傳播方向〞，由于在z=0和z=l處存在導體壁，電磁波將在期間來回反射形成駐波，所以在空腔內(nèi)不可能有波的傳播。3.3三維問題〔直角坐標系〕先討論式第一個公式：--(7)在時間步，對節(jié)點的離散公式為：--(8)上式中的第二項用平均值來替代是因為離散方程中電場的時間取樣是整數(shù)n，磁場的時間取樣是n+1/2，所以只能取n及n+1時電場的平均值。實際也證明這個平均值使FDTD算法具有數(shù)值穩(wěn)定性。整理后，將作為未知數(shù)，其余作為迭代計算的數(shù)式中：--〔9〕同理，式〔1-3〕中其它兩個公式的離散形式為式中：--〔10〕式中：--〔11〕以上三式是電場的時間推進計算公式。同樣，討論式〔3-4〕中第一個公式，設觀察點為的節(jié)點，即在時刻，對節(jié)點的離散公式為：式中，----〔12〕同理，式〔1-4〕中其它兩個公式的離散形式為式中，--〔13〕式中，--〔14〕以上三式是磁場的時間推進計算公式。時域推進計算框圖〔交叉半步逐步推進〕假設時空間各節(jié)點處的電場值假設時空間各節(jié)點處的電場值〔賦初值〕計算時空間各節(jié)點處的磁場值，式〔8〕—〔11〕計算時空間各節(jié)點處的電場值，式〔8〕—〔10〕在編程中，為了使電場和磁場有相同的數(shù)量級〔為減小誤差〕，可對H或E進行“歸一化〞處理，即：用取代，用取代，式中是自由空間波阻抗。計算結果再分別除以和乘以即可?？梢钥闯觯@種離散方法電場和磁場在時間順序上交替抽樣，抽樣間隔相差半個時間步長，使麥克斯韋方程離散后成為顯示差分方程，從而可以在時間上迭代求解，不需矩陣求逆。給定初值后，可以逐步推進，求得以后各個時刻點的空間電磁分布。這是FDTD法的最大特點。3.4FDTD法的穩(wěn)定性及數(shù)值色散問題在FDTD中，時間增量和空間增量、、之間不是相互獨立的，它們的取值必須滿足一定的關系，以防止數(shù)值結果的不穩(wěn)定——表現(xiàn)為隨著時間步數(shù)的增加，計算結果發(fā)散。造成解不穩(wěn)定的因素有多種：計算機在計算過程中，原始數(shù)據(jù)可能有誤差，如系數(shù)陣建立過程中產(chǎn)生的誤差，而每次運算由于只能保存有限位數(shù)而又產(chǎn)生誤差，誤差的積累有可能淹沒真正解，使計算結果不可靠，即不穩(wěn)定；計算方法不適宜；離散間隔不當?shù)取榱舜_定數(shù)值解穩(wěn)定的條件，有許多推導方式，結論相同.1、時間步長穩(wěn)定性要求〔推導略〕一般情況下〔15〕T為波動周期。2、時間步長與空間步長的關系三維〔16〕在非均勻區(qū)域，v取最大值。真空中v=c(光速)。假設是立方體Yee元胞，，那么〔17〕。假設是正方形Yee元胞〔二維〕，，那么〔18〕。假設是線段等分Yee元胞〔一維〕，那么〔19〕假設電磁波所在空間的介質(zhì)特性與頻率有關，則電磁波的傳播速度也將是頻率的函數(shù)，這種現(xiàn)象稱色散。而FDTD方程是原Maxwell方程的一種近似，所以當計算機對電磁波在空間的傳播進行模擬時，在非色散空問中也會出現(xiàn)色散現(xiàn)象，且電磁波的相速度隨波長、傳播方向及變量離散化的情況發(fā)生變化，這種非物理性的色散稱為數(shù)值色散。數(shù)值色散會導致脈沖波形的破壞，出現(xiàn)人為的各向異性及虛假的折射現(xiàn)象。數(shù)值色散是由于近似差商替代連續(xù)微商引起的，這種影響可以通過減小離散化過程所取空間和時間步長而無限減小，但計算空間的總網(wǎng)格數(shù)目的增加也會相應增加對計算機存儲空間和計算時間的要求，所以要根據(jù)實際條件來選擇適宜的步長。而數(shù)值色散在FDTD法中是不可防止的。為了減小數(shù)值色散，除了滿足式〔16〕外，可取比式〔15〕更高的要求〔20〕并令電磁波沿網(wǎng)格的對角線方向傳播，這時有，可以得到理想的色散關系。3.5建模及原程序%%矩形腔寬度為25mm,高度為12.5mm,諧振腔長度60mm，介質(zhì)為空氣%%計算得諧振腔的諧振頻率為13.42Ghzclear;clc;c=3e8;muz=4.0*pi*1.0e-7;epsz=1.0/(c^2*muz);%每小格長,三個方向的網(wǎng)格長均設為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;%總步數(shù)nmax=[800,1000，3000];,%%選擇f=1;nmax=nmax(f);%脈沖的系數(shù)rtau=0.3/f0;%決定脈寬的時間長度tau=rtau/dt;%歸dt化ndelay=tau*3;%時延%介質(zhì)的參數(shù)%選擇腔內(nèi)任意一點進行驗證i=10;j=5;k=10;%介質(zhì)的參數(shù)epsr=1.0;sig=0.0;mur=1.0;sigm=0.0;%迭代系數(shù)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結果及分析矩形波導諧振腔是由兩端短路的一段金屬波導構成，矩形波導參數(shù)----寬度a=25mm，高度b=12.5mm，長度l=60mm。由微波理論可知，當b<a<l時，TE101模的諧振波長最長，是矩形波導諧振腔的主模。TE波〔即橫電波〕是在波的傳播方向上有磁場分量，但沒有電場分量，即電場垂直于傳播方向。TE101模在矩形腔的3個方向都不傳輸能量，呈駐波分布，并且由邊界條件決定：電場分量只有Ez分量，在腔體中央最強；磁場有Hx和Hy兩個分量，在腔壁附近最強，腔體中央為0。由于矩形波導諧振腔有著廣泛的應用，并且大多工作在主模狀態(tài)，這里選擇微分高斯脈沖以鼓勵起模?！?〕決定網(wǎng)格單元的尺寸和時間步長取Δx=Δy=Δz=Δ=1.25mm，則空間網(wǎng)格數(shù)為50×25×120，取Δt=Δ/(2c)=2.0833ps?！?〕設置邊界條件對于矩形波導諧振腔，腔體的6個面都是金屬，并且均為理想導體，即腔內(nèi)導體邊界上的所有切向電場分量為0，所有法向磁場分量為0。矩形諧振腔的諧振波長為：QUOTE，則TE101模的諧振波長為=22.3606mm，主模諧振頻率為=13.42Ghz〔3〕設置鼓勵源高斯脈沖的表達式為:QUOTEQUOTE為了在諧振腔中鼓勵起模，并且抑制其他高次模，選擇線源脈沖，使之在腔內(nèi)xy平面中心處沿y軸方向分布，并選擇適宜的和τ值。經(jīng)過反復試驗，取=138.462ps，τ=46.154ps。圖4-1分析時長為〔dt*nmax=2.0833ps*800〕注意對應時，在諧振腔體中上圖顯示了在諧振腔內(nèi)部由高斯脈沖源鼓勵起由諧振腔自身參數(shù)選擇出的諧振頻率分量，〔12.5，6.25，12.5〕mm處各個長分量所包含的頻率分量與理論值比照，發(fā)現(xiàn)在諧振頻率13.2Ghz處出產(chǎn)生了最大響應。波形圖中的其他波峰是由于高斯鼓勵脈沖元不是理想脈沖所引起的，對本次驗證試驗沒有太大影響。因為不是絕對理想情況，與計算所得的13.42數(shù)據(jù)誤差相差不大，即得到驗證。4〕時間和空間離散度以及采樣點數(shù)對分析結果的影響時間離散度不但應滿足數(shù)值色散的最低要求，而且要與空間補償滿足穩(wěn)定關系，在這次試驗中，根據(jù)Yee元胞模型，在三個坐標方向均采用等步長dx=1.25mm，根據(jù)3-12,3-15計算出滿足條件的是