基于二維有限差分法對非對稱微帶線TM模的計算分析

基于二維有限差分法對非對稱屏蔽帶狀線TM模的計算分析摘要：利用二維有限差分法對非對稱屏蔽帶狀線的TM模進行計算。通過調(diào)整網(wǎng)格數(shù)目觀察計算結(jié)果的精度，可以設(shè)定出最優(yōu)的網(wǎng)格數(shù)目，然后改變微帶尺寸，觀察最低兩個TM模對應(yīng)的截止頻率的變化規(guī)律，最后給出對應(yīng)的橫截面電場、磁場分布圖。關(guān)鍵詞：有限差分法，非對稱屏蔽帶狀線，TM模TheCalculationandAnalysisofUnsymmetricalShieldedStrip-LineBasedon2-DimensionFDMAbstract:TheTMmodeofUnsymmetricalShieldedStrip-Lineiscalculatedusing2-dimensionFDM.Byadjustingthenumberofmeshgrid,accordingtotheprecisionofthecalculatedresult,theoptimummeshgridnumberisobtained.Then,bychangingthesizeoftheUnsymmetricalShieldedStrip-Line,thecharacterofthelowesttwoTMmodeisanalyzed.Finally,thechartoftheelectricfieldandmagneticfieldonthecross-sectionisdrawn.Keyword：FDM,UnsymmetricalShieldedStrip-Line,TMMode赫姆赫茲方程及邊界條件圖1屏蔽帶狀線示意圖圖1所示的屏蔽帶狀線內(nèi)傳播的TM波的Ez分量，滿足赫姆赫茲方程(1)及邊界條件：〔當X、Y位于內(nèi)外邊界上時〕。求出后就可畫出最低兩個TM模的色散特性曲線。色散方程為〔7〕而求出后可根據(jù)以下公式求橫向電磁場分量：(2)(3)(4)(5)2.離散化場域、建立差分模型并構(gòu)造稀疏矩陣K圖2五點差分格式示意圖參考文獻1中的方法，采用5點差分格式〔圖2〕，將場域進行離散化處理，對圖1中的帶狀線進行網(wǎng)格剖分。此時電場縱向分量Ez在內(nèi)、外邊界上滿足Ez＝0的邊界條件。以A〔i,j〕代表縱向電場Ez在網(wǎng)格各點的的值。對〔1〕式進行離散化，選取步長hx＝hy＝h。那么：〔6〕帶狀線上的點的差分方程：緊靠帶狀線的非金屬點滿足的差分方程為帶狀線同一層其他非金屬點滿足的差分方程：帶狀線正上方的點的差分方程:帶狀線正下方的點的差分方程:緊靠左側(cè)波導(dǎo)壁的一列：緊靠右側(cè)波導(dǎo)壁的一列緊靠底部波導(dǎo)壁的一行緊靠頂部波導(dǎo)壁的一行將所有點的差分方程寫出來，并把項全部移到等號右邊，得其中，，由此可以計算出矩陣的特征值〔〕，每個都對應(yīng)一個特征頻率，即一個TM模，最低的兩個非零的特征頻率對應(yīng)最低的兩個TM模。特征值對應(yīng)的特征向量就是該模式的縱向電場Ez。根據(jù)式〔2〕－〔5〕便可求得橫截面的電場、磁場分布。3.算法流程圖根據(jù)第2節(jié)建立的模型，可以設(shè)計初如下算法，見圖4。啟啟動給出差分方程組的系數(shù)矩陣A求矩陣A的特征值及所對應(yīng)的特征矢量選出最小的兩個非零特征值及所對應(yīng)的特征矢量將特征矢量轉(zhuǎn)化為場域上按節(jié)點排列的Ez矩陣求出橫截面電場和磁場所對應(yīng)的矩陣并畫出電磁場圖求出幾個最低模式的截止波數(shù)截止波長，畫出所對應(yīng)的色散特性曲線結(jié)束圖4算法流程圖表1.TM模截至頻率隨網(wǎng)格數(shù)目劃分的變化表格網(wǎng)格1216202428323640第1個特征頻率〔*10GHz〕2.45572.46542.47032.47322.47522.46472.46742.4695第2個特征頻率〔*10GHz〕3.49003.51723.53083.53883.54393.52203.52813.53284.計算結(jié)果分析：4.1．網(wǎng)格粗細對結(jié)果的影響最低兩個TM模的截止頻率的計算值隨網(wǎng)格粗細的變化情況如表1所示。從圖中可以看出，截止波長計算值隨網(wǎng)格數(shù)的增加而趨于穩(wěn)定，網(wǎng)格越細越接近真實值，但網(wǎng)格越細計算量越大，耗時越長；由上分析可以得出取網(wǎng)格數(shù)為20來處理問題，已到達比擬高的精度。由此，確定網(wǎng)格數(shù)目為20，對屏蔽帶狀線作進一步分析表2所編程序計算結(jié)果與CST仿真結(jié)果的比照Freq1〔GHz〕Freq2〔GHz〕程序結(jié)果24.70335.308CST計算結(jié)果25.26936.3954.2計算結(jié)果與CST仿真的結(jié)果進行比擬為了驗證程序是否正確，取a/b=0.5，計算前2個截止頻率與CST仿真結(jié)果比照，如表2所示可見，程序計算結(jié)果與CST仿真結(jié)果比擬一致。4.3．色散曲線取a/b＝0.5時，可以得到最低兩個TM模色散曲線，如圖5所示。而后通過改變a/b的值，可以得到最低兩個TM模對應(yīng)的特征頻率隨a/b的變化的情況，，見圖6。圖5色散曲線圖6最低兩個TM模對應(yīng)的特征頻率隨a/b的變化率曲線圖4.4.橫截面上電場，磁場分布圖通過以上分析，根據(jù)〔1〕、〔2〕－〔5〕式最終可以繪制出截面上電場，磁場分布圖，見圖7。7〔a〕a/b＝0.3，最低TM模與次低TM模的電、磁場7〔b〕a/b＝0.5，最低TM模與次低TM模的電、磁場7〔c〕a/b＝0.7，最低TM模與次低TM模的電、磁場結(jié)論由赫姆赫茲方程及邊界條件，利用二維有限差分法對屏蔽帶狀線TM模計算。在計算出縱向電場分量后，進而求出電場、磁場的其它分量。然后對屏蔽帶狀線的特性進行了分析，通過觀察數(shù)值計算結(jié)果，發(fā)現(xiàn)網(wǎng)格數(shù)目在20時，比擬恰當。并由此給出橫截面電場、磁場的分布圖。參考文獻[1]王秉中.計算電磁學(xué).科學(xué)出版社,2002附錄：程序代碼主程序%非對稱微帶線clearallcloseall%clfclc%formatshortglobalNrhomueepsiloncbahtime=cputime;N=20;H=N;%N是寬度，H是高度c=3e8;rho=0.7;b=5e-3;h=2*b/N;%波導(dǎo)邊長分為N段，網(wǎng)格步長為ha=rho*b;M=round((b-a)/h);rho1=0.5;%高度比例因子high=b+rho1*b;mue=pi*4e-7;epsilon=1/(pi*36e9);high=round(high/h);%這個h1是帶狀線的高度%H=b+mesh*h;%嘗試把微帶線上調(diào)5個網(wǎng)格%構(gòu)造I矩陣I=eye(N-1);%構(gòu)造D矩陣.這是除微帶線以外的A〔i，j〕的稀疏矩陣%矩陣的構(gòu)造主要是對主對角線，以及上、下一層的數(shù)據(jù)的設(shè)定,D矩陣是除了微帶線這一層外，其他的對角矩陣塊。D(1:N-1,1:N-1)=0;fori=1:N-1D(i,i)=-4;endfori=1:N-2D(i,i+1)=1;D(i+1,i)=1;end%構(gòu)造stripline矩陣.?當N/2-M-1>i時，會有微帶接近邊界的情況，此時stripline矩陣就會有所變化...?%此時的stripline矩陣是微帶線這一層的矩陣，不管怎么變化，主對角線的值是保持－4不變化的。fori=1:N-1stripline(i,i)=-4;endfori=2:N-2stripline(i,i-1)=1;stripline(i,i+1)=1;endi=N/2-M-1stripline(i,i-1)=1;stripline(i,i+1)=0;i=N/2+M+1stripline(i,i-1)=0;stripline(i,i+1)=1;fori=(N/2-M):(N/2+M)stripline(i,i-1)=0;stripline(i,i+1)=0;endstripline(1,2)=1;stripline(N-1,N-2)=1;%up_down矩陣這是微帶線上一層A(i,j-1)，下一層的A(i,j+1)的網(wǎng)格的系數(shù)。是否為零？fori=1:N-1ifi>=(N/2-M)&&i<=(N/2+M)up_down(i,i)=0;elseup_down(i,i)=1;endend%%%。。。。。。。。。。。K矩陣。。。。。。。。。。。。。。。。。。。。。。。K(1:(N-1)*(N-1),1:(N-1)*(N-1))=0;fori=2:N-2ifi==high-1K(((i-1)*(N-1)+1):(i*(N-1)),((i-2)*(N-1)+1):((i-1)*(N-1)))=I;K(((i-1)*(N-1)+1):(i*(N-1)),((i-1)*(N-1)+1):(i*(N-1)))=D;K(((i-1)*(N-1)+1):(i*(N-1)),(i*(N-1)+1):((i+1)*(N-1)))=up_down;elseifi==highK(((i-1)*(N-1)+1):(i*(N-1)),((i-2)*(N-1)+1):((i-1)*(N-1)))=up_down;K(((i-1)*(N-1)+1):(i*(N-1)),((i-1)*(N-1)+1):(i*(N-1)))=stripline;K(((i-1)*(N-1)+1):(i*(N-1)),(i*(N-1)+1):((i+1)*(N-1)))=up_down;elseifi==high+1K(((i-1)*(N-1)+1):(i*(N-1)),((i-2)*(N-1)+1):((i-1)*(N-1)))=up_down;K(((i-1)*(N-1)+1):(i*(N-1)),((i-1)*(N-1)+1):(i*(N-1)))=D;K(((i-1)*(N-1)+1):(i*(N-1)),(i*(N-1)+1):((i+1)*(N-1)))=I;elseK(((i-1)*(N-1)+1):(i*(N-1)),((i-2)*(N-1)+1):((i-1)*(N-1)))=I;K(((i-1)*(N-1)+1):(i*(N-1)),((i-1)*(N-1)+1):(i*(N-1)))=D;K(((i-1)*(N-1)+1):(i*(N-1)),(i*(N-1)+1):((i+1)*(N-1)))=I;endendK(1:N-1,1:N-1)=D;K(1:N-1,N:2*(N-1))=I;K((N-2)*(N-1)+1:(N-1)^2,(N-3)*(N-1)+1:(N-2)*(N-1))=I;K((N-2)*(N-1)+1:(N-1)^2,(N-2)*(N-1)+1:(N-1)^2)=D;%..................................................................%[phi,value]=eig(K);%求出特征向量矩陣phi，特征值陣eigv[phi,value]=eigs(K,10,'SM');eig_freq=c*sqrt(-value)/(2*pi*h);%根據(jù)特征值，求出特征頻率，[freq1,position1,freq2,position2]=min_eigvalue(eig_freq);phi1=phi(:,position1);phi2=phi(:,position2);f=freq1:0.5e9:50e9;beta1=2*pi*sqrt(mue*epsilon*(f.^2-freq1.^2));beta2=2*pi*sqrt(mue*epsilon*(f.^2-freq2.^2));figure(1)plot(f,beta1,'*',f,beta2,'--')axisnormaltitle('freq-beta曲線');xlabel('頻率freq');ylabel('相移常數(shù)beta');legend('最低TM模f-beta曲線','次低模f-beta曲線')%...................求出縱向場之后，再求出橫向分量.........................kc=2*pi*freq1/c;w=2*pi*freq1;Z0=sqrt(mue/epsilon);%TM模的Hy分量%phi3=phi(:,121);%phi1=phi3';Ez(1:N-1,1:H-1)=0;forn=1:H-1Ez(n,1:N-1)=phi1((n-1)*(N-1)+1:n*(N-1));end%forn=1:H-1%Ez(1:N-1,n)=phi1((n-1)*(N-1)+1:n*(N-1))';%end%...................求出橫向分量............................j=j;fork=1:H-1fori=2:N-2Hy(i,k)=(-j*w*epsilon/(2*h*kc^2))*(Ez(i+1,k)-Ez(i-1,k));endHy(1,k)=(-j*w*epsilon/(2*h*kc^2))*Ez(2,k);Hy(N-1,k)=(j*w*epsilon/(2*h*kc^2))*Ez(N-2,k);end%TM模的Hx分量fori=1:N-1fork=2:H-2Hx(i,k)=(j*w*epsilon/(2*h*kc^2))*(Ez(i,k+1)-Ez(i,k-1));endHx(i,1)=(j*w*epsilon/(2*h*kc^2))*Ez(i,2);Hx(i,H-1)=-(j*w*epsilon/(2*h*kc^2))*Ez(i,N-2);endEx=Z0*Hy;Ey=-Z0*Hx;figure(2)subplot(2,2,1)%quiver((Ey),(Ex));Efield1=quiver(imag(Ey),imag(Ex));%第1個TM模的橫截面電場分布圖title('TM模的橫截面電場分布圖')subplot(2,2,2)%quiver((Hy),(Hx));Hfield1=quiver(imag(Hy),imag(Hx));%第1個TM模的橫截面磁場分布圖title('TM模的橫截面磁場分布圖')Ez(1:N-1,1:H-1)=0;forn=1:H-1Ez(n,1:N-1)=phi2((n-1)*(N-1)+1:n*(N-1));end%forn=1:H-1%Ez(1:N-1,n)=phi1((n-1)*(N-1)+1:n*(N-1))';%end%...................求出橫向分量............................j=j;fork=1:H-1fori=2:N-2Hy(i,k)=(-j*w*epsilon/(2*h*kc^2))*(Ez(i+1,k)-Ez(i-1,k));endHy(1,k)=(-j*w*epsilon/(2*h*kc^2))*Ez(2,k);Hy(N-1,k)=(j*w*epsilon/(2*h*kc^2))*Ez(N-2,k);end%TM模的Hx分量fori=1:N-1fork=2:H-2Hx(i,k)=(j*w*epsilon/(2*h*kc^2))*(Ez(i,k+1)-Ez(i,k-1));endHx(i,1)=(j*