1、電磁仿真綜合實踐Numerical Simulation of Electromagnetic Field作者 張泳 Author Yong Zhang指導老師 王秉中Instructor Bing-Zhong Wang單位 電子科大物理電子學院99777Class 99777 CPE of UESTC設計時間 2002年5月8日6月13日【摘要】仿真是通過建立系統模型，并且對模型進行實際研究，以代替對于實際系統的研究。Matlab 5.3是一種功能強大工程數學軟件。本文闡明了怎樣利用Matlab 5.3實現一維以及二維電磁仿真實例的建模分析、編程調試、結果分析。【Abstract】Simul

2、ation is an available means in studying real system by Model-building.This paper building a numerical simulation of electromagnetic field model based on the Matlab 5.3,the design is discussed in detail.【關鍵詞】電磁仿真，計算電磁學，電磁場【Key Words】 Electromagnetic Simulation,Computational Electromagnetics,Electroma

3、gnetic Field【設計目的】1 課程性質：屬于非計算機專業計算機基礎教育第三層次計算機綜合應用層次的課程2 任務： 21 加強學生的計算機綜合應用能力、尤其是運用計算機分析和解決專業問題的能力培養，使學生對獨立進行科學研究有初步實踐。 22 初步掌握一種純數值電磁仿真方法時域有限差分法。 23 初步學會綜合應用一種程序設計語言進行科學與工程計算。 24 增強科技報告的寫作能力，學會相關軟件使用。正文第一部分：電磁仿真的意義現代科學研究的基本模式：“科學實驗，理論分析，數值計算”三位一體。計算電磁學：以電磁場理論為基礎，以計算機及計算機技術為工具和手段，運用計算科學所提供的各種方法，解決

4、復雜電磁問題的一門應用科學。電磁場工程電磁場理論 電磁規律 實驗驗證 數學方程計算電磁學 計算數據 建模與仿真 運算手段和方法 優化與設計【電場問題求解方法】解析法電磁場問題的傳統解法。強調電磁分析和數學分析通常給出緊湊高效的計算程序。程序的最終擁護只具有很少的彈性，只能改變很少的參數，結構的主要特性已被編進程序，適用于專用程序開發。純數值法直接以數值的、程序的形式描述電磁場問題。普適性強，用戶彈性大。特定問題的邊界條件、電氣結構、激勵等特性可以不編入基本程序，而由用戶輸入，更好的情況是通過圖形界面輸入。用戶不必具備高深的電磁理論、數學及數值技術方面的知識。適用于通用商業程序開發。受到硬件限制

5、大。 FDTD TLM FEM 惠更斯原理 變分原理麥克斯韋方程數學描述 及 離散化處理結果計算計算程序數值模型 邊界條件材料特性 后處理 用戶界面 第二部分 時域有限差分法基本原理【差分原理】前項差分 誤差 O(h)后項差分 誤差 O(h)中心差分 誤差 O(h*h) 用中心差分代替對空間、時間坐標的微分 【Yee的差分網格】每個磁場分量被四個電場分量環繞，每個電場分量被四個磁場分量環繞。E和H的空間取樣位置相差半個空間步長，E和H時間取樣位置相差半個時間步長。Yee的差分算法考慮空間一個無緣區域,其媒質的參數不隨時間變化且各項同性,則麥克斯韋旋度方程在直角坐標系中,寫成分量式:= (1)

6、(2) (3) (4) (5) (6)【蛙跳格式】H0 H1/2 H3/2E0 E1 E2【解的穩定條件】：； 【數值色散】 在時域有限差分網格中，數值波模的傳播速度將隨頻率改變，即有色散。這種色散由數值網格引起，而非物理上客觀存在。為了減小色散，實用中通常取空間步長滿足 。第三部分 均勻平行板傳輸線傳輸特性仿真 FDTD基本過程【平行板傳輸線中的主模TEM模】根據平行板傳輸線中TEM波的特性對MAXWELL旋度方程簡化得： 其中 若取時域波形為高斯波形， ，如下圖所示。 【一維差分格式】E【網格劃分】Ex Hy Ex Hy Ex Hy Ex Hy Ex k=0 k=1/2 k=1 k=3/2

7、 k=2 k=Nz-1/2 k=Nz T=0.5ns，則=1GHz，d=0.18m,L=6m, =0.3m, =0.015m, 400, 【時間步進過程】在初始（n=0）時令所有場兩為零。將高斯脈沖激勵加在輸入口 由FDTD方程算出n+1/2時刻所有場點的磁場 由FDTD方程算出n+1時刻所有場點的電場 令金屬邊界上的切向電場分量為零，并在截斷邊界上使用截斷邊界條件 記錄并輸出兩個觀測面處的電壓值 n n+1重復上一步驟，直到（n=nt）脈沖波全部傳出計算區域。 【截斷邊界條件】 終端短路： 終端匹配：（吸收邊界條件） 【一維仿真實例】終端匹配：入射波大約經過1000步后到達邊界，由于采用吸收

8、邊界條件，高斯波幾乎被完全吸收，通過圖可以觀測到吸收后只有振幅小于的波形存在，這在誤差的允許范圍之內。第一條、第二條分別為第一、二觀測面的電壓。時間步長sZ方向網格數（0015m）Z方向網格數（0015m）電壓強度V電場強度V/M電場強度V/M第二條第一條 入射波 被吸收后的波形 觀測面電壓終端短路：入射波到達終端后由于截斷邊界條件采取的是短路，所以波形反轉，向入射口傳輸，由于Nt取得足夠大（Nt=2000），波形會到達入射口，所以應該在波形到達入射口之前，激勵脈沖已完全進入入射口后，將入射口邊界條件切換為吸收邊界條件。第一條第二條分別為第一、二觀測面的電壓。在橫軸下方的是反射電壓，上方為入射

9、電壓。第二條電場強度V/M時間步長sZ方向網格數（0015m）Z方向網格數（0015m）第一條第一條電壓強度V電場強度V/M 被反射的波形 被吸收后的波形 觀測面電壓第四部分 帶擋板的平行板傳輸線傳輸特性仿真X L/2 d s=d/3Z L【二維FDTD方程】由于不連續性的出現，沿X方向場的均勻性遭到破壞，沿Y方向場的均勻性仍然保持，此時平行板傳輸線中的場分量有Ex,Ez,Hy.【參數設置】：數組EX(NX,NZ+1),EZ(NX+1,NZ),HY(NX,NZ),L/2應選足夠長，使得在擋板反射波返回入射口之前，激勵脈沖已完全進入入射口，并且入射口邊界條件已切換為吸收邊界條件。【觀測面電壓波形

10、】【入射、反射、透射電壓頻譜】【網絡的S參數】現在我們再記錄下總電壓、就可以分離出入射電壓、反射、透射電壓。=- = 終端短路的反射系數：終端匹配的S參數： 【終端條件】終端匹配終端短路=0【觀測面電壓波形】【入射、反射、透射電壓頻譜】【二維仿真實例】 終端匹配：可以通過對Ex波形的觀測知道，波形到達擋板時產生尖端效應，高斯波有一部分被“擠”過擋板，有一部分被擋板反射。此時應該將入射口和終端都設為吸收邊界條件。 如下圖所示：Z方向網格數(0.015m)X方向網格數(0.015m)電場強度Ex V/MZ方向網格數(0.015m)X方向網格數(0.015m)電場強度Ex V/M X方向網格數(0.

11、015m) 入射波 到達擋板 電場強度Ex V/MZ方向網格數(0.015m)Z方向網格數(0.015m)X方向網格數(0.015m)X方向網格數(0.015m)電場強度Ex V/M 通過擋板 兩端都吸收電壓幅度(V)時間步長（s）第一個觀測面反射電壓第二個觀測面透射電壓第一個觀測面入射電壓 參考面電壓 S參數頻率f (0.1GHz)終端匹配時S參數終端短路：由于擋板與短路板之間構成一種諧振腔結構，時域脈沖波在其間來回反射，入射口仍然是吸收邊界條件，波完全消失的時間會比終端匹配情況要長。Nt足夠大（Nt=20000）如下圖所示：電場強度ExV/MX方向網格數(0.015m)z方向網格數(0.0

12、15m) 終端短路、入射口吸收，Ex被終端反射的波形電壓幅度V 頻率f (0.1GHz) 終端短路時第一參考面電壓 反射系數頻率f (0.1GHz) 終端短路時反射系數 第五部分 課程總結與討論短短幾周的電磁仿真綜合實踐很快就結束了，這幾周緊張的學習給我留下了深刻的印象。我們已經學習了電磁場與電磁波、微波技術基礎等課程，對電磁理論有一定的認識和了解，但那只是完全理論性質的。通過這門課程的學習，我對電磁場等相關知識有了進一步的了解，有了一定感性的認識，感覺電磁場也不是那么和天書一樣的難以理解。通過用計算機模擬將電磁現象展示在我們面前，使得我不但加深了對電磁理論的理解，同時也鍛煉了我們計算機應用方

13、面的能力，對科學研究有了一定的了解。我有如下幾點體會：首先，這次實踐使得我對計算機的應用有了進一步的提高。以往我對計算機的認識只是處于一個很簡單的層次上，對其應用、原理也只是略知皮毛。但通過這段時間的學習，自己用計算機以及相關軟件完成了這次實踐活動。雖然在開始的時候對軟件以及原理理解不夠，但是通過老師的多次講解，同學之間的互相探討，漸漸地自己掌握了一定的方法。看著自己編出的程序調試出令人滿意的結果，還是相當喜悅的。其次，這次的實踐活動讓我加深了對專業知識的了解，對電磁理論也不象以前那樣見到就躲了。當我用自己編的程序捕捉到了它們，對它們的產生、傳播、反射、透射等現象有了更為直接的認識和了解。自己

14、本來對電磁理論有一定的興趣，苦于沒有很好的方法幫助自己更好的理解和學習，或許通過這次實踐活動，自己以后能夠找到適合自己的學習電磁理論的方法。第三，通過對原理、程序的不斷學習、操作，使得我認識到做任何事情都要認真，持之以恒。王老師科學嚴謹的作風給我留下了深刻的印象，不論是做人還是做科研，王老師都是那么的一絲不茍。這一點讓我也不敢有所懈怠。同時我也感覺到自己能力的欠缺，分析問題、解決問題的能力還有待提高。相信在以后我會通過自己的不斷努力，提高自己的能力，學習和工作上更加細致認真。第四，這次實踐活動，我感覺是我上大學以來收獲最大的一門課程。王老師講課不是傳統的灌輸式的教育，學生是真正意義上學習的主體

15、。課前需要我們對理論方面的知識多加復習和學習，上課的時候需要我們積極的開動腦筋，如何把理論轉化為實際可見，課后也需要我們不斷的思考，如何才能使得我們做的軟件更加趨于完善。在掌握計算電磁學有關理論知識的同時，更加重要的是我們學到了方法，這一點我覺得才是最為重要的。下面說一下我對這門課程的一些建議吧。我覺得要消除一些學生對這門課程的緊張情緒（開始我就有這樣的感覺），適當的增加課時數，畢竟我們目前涉及的真的只是皮毛而已。增加課時數，可以讓大家對理論的理解加深，同時也會有更多的同學以后投入到電磁場領域方面的學習和研究。還有上機時間最好不要靠得太近了，畢竟學習有一個消化吸收的過程。這次的科研實踐活動是令

16、人愉快的。它讓我體會到了快樂，當然也有抓耳撓腮的痛苦；它讓我學習到了更多新的知識，新的方法；它讓我看到了自己的不足之處，在以后的學習過程中，我一定會加以改正的。最后我還要特別感謝王老師，以及研究生師姐在這次實踐活動中給予我的幫助。正是他們一絲不茍的工作態度和嚴謹的作風深深的教育了我，以后我也會盡我最大的努力去做好每一件事情。第六部分 附錄一、參考文獻1、計算電磁學（節選），王秉中，電子科技大學，2000年版2、電磁仿真綜合實踐，王秉中，電子科技大學3、MATLAB基礎與應用簡明教程，張平 ，北京航空航天大學出版社，2001年1月第一版二、相關程序1、一維終端匹配%定義變量t=0.5E-9;c=

17、3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*p;y0=120*pi;dz=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=1000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cmur=(dt*c-dz)/(dt*c+dz);ex(1:nz+1)=0;hy(1:nz)=0;vref1(1:nt)=0;vref2(1:nt)=0;%程序for n=1:nt ex(1)=exp(-(n*dt-t0)2/t2); for i=1:nz

18、 hy(i)=hy(i)+cm*(ex(i)-ex(i+1); end exm=ex(nz); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); end ex(nz+1)=exm+cmur*(ex(nz)-ex(nz+1);%一維匹配終端吸收條件 vref1(n)=ex(k1)*d; vref2(n)=ex(k2)*d; %plot(1:nz+1,ex);pause(0.001);%電場強度監控 x=1:nt; plot(x,vref1(x),'r',x,vref2(x),'g');pause(0.001);%觀測面電壓圖 %g

19、ridend%一維終端匹配2、一維終端短路%定義變量y0=120*pi;dz=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=2000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cmur=(dt*c-dz)/(dt*c+dz);ex(1:nz+1)=0;hy(1:nz)=0;vref1(1:nt)=0;vref2(1:nt)=0;for n=1:nt for i=1:nz hy(i)=hy(i)+cm*(ex(i)-ex(i+1); end if n<=800 ex(1)=exp(-(n

20、*dt-t0)2/t2); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); endelse exm=ex(2); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); end ex(nz+1)=0; %終端短路邊界條件 ex(1)=exm+cmur*(ex(2)-ex(1); %高斯波完全進入后，入口改為吸收邊界條件end vref1(n)=ex(k1)*d; vref2(n)=ex(k2)*d; x=1:nt;plot(x,vref1(x),'r',x,vref2(x),'g');%兩個觀測面的

21、電壓波形 pause(0.01);% plot(1:nz+1,ex);pause(0.001);%實時監控Ex波形 %gridend%一維終端短路3、二維終端匹配%定義變量t=0.5E-9;c=3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*pi;y0=120*pi;dz=lammin/20;dx=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=1000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cf=dt/(epslon*dx)

22、;cmur=(dt*c-dz)/(dt*c+dz);nx=12;ez(1:nx+1,1:nz)=0;ex(1:nx,1:nz+1)=0;hy(1:nx,1:nz)=0;ver1in(1:nt)=0;ver1re(1:nt)=0;ver2tr(1:nt)=0;exb(1:nx,1)=0;exe(1:nx,1)=0;nf=10;df=1e8;jj=sqrt(-1);gver1in(1:nf)=0;gver1re(1:nf)=0;gver2tr(1:nf)=0;s11(1:nf)=0;s21(1:nf)=0;for n=1:nt %時間循環 for i=1:nx %二維FDTD過程 for k=1:

23、nz hy(i,k)=hy(i,k)+cm*(ez(i+1,k)-ez(i,k)*dz/dx+ex(i,k)-ex(i,k+1); end end if n<=300 ex(1:12,1)=exp(-(n*dt-t0)2/t2); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; %設置擋板 for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy

24、(i,k)-hy(i-1,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; else %高斯波完全進入，入射端改為吸收邊界條件 exb(1:nx,1)=ex(1:nx,2);exe(1:nx,1)=ex(1:nx,nz); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy(i,k)-hy(i-1,k); end endex

25、(1:nx,1)=exb(1:nx,1)+cmur*(ex(1:nx,2)-ex(1:nx,1);%入口邊界條件ex(1:nx,nz+1)=exe(1:nx,1)+cmur*(ex(1:nx,nz)-ex(1:nx,nz+1);%吸收邊界條件 ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0;endfor i=1:nx %分離出入射電壓、反射電壓、透射電壓 ver2tr(n)=ver2tr(n)+ex(i,301)*d/12; end if n<=500 for i=1:nx ver1in(n)=ver1in(n)+

26、ex(i,101)*d/12; end else for i=1:nx ver1re(n)=ver1re(n)+ex(i,101)*d/12; end end %x=1:nt%plot(x,ver1in(x),'r',x,ver1re(x),'b',x,ver2tr(x),'m');%各個電壓波形surf(ex); %監控Ex波形pause(0.01)endfor f=1:nf for n=1:nt %傅氏變換 gver1in(f)=gver1in(f)+dt*ver1in(n)*exp(-2*pi*f*df*n*dt*jj); gver1re(

27、f)=gver1re(f)+dt*ver1re(n)*exp(-2*pi*f*df*n*dt*jj); gver2tr(f)=gver2tr(f)+dt*ver2tr(n)*exp(-2*pi*f*df*n*dt*jj); endendfor f=1:nf %S參數 s11(f)=gver1re(f)/gver1in(f); s21(f)=gver2tr(f)/gver1in(f);endfor i=1:nf s1(i)=abs(s11(i); s2(i)=abs(s21(i);end%x=1:nf;%plot(x,s1,'g',x,s2,'r');%S參數%g

28、rid %二維終端匹配4、二維終端短路%定義變量t=0.5E-9;c=3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*pi; y0=120*pi;dz=lammin/20;dx=lammin/20;dt=dz/(2*c); t0=3*t;l=6;d=0.18;nz=400; nt=20000;k1=101;k2=301; cm=dt/(mu*dz);ce=dt/(epslon*dz); cf=dt/(epslon*dx);cmur=(dt*c-dz)/(dt*c+dz);nx=12; ez(1:nx+1,1:nz)=0;ex(1:

29、nx,1:nz+1)=0;hy(1:nx,1:nz)=0; vref1(1:nt)=0;vref2(1:nt)=0;ver1in(1:nt)=0; ver1re(1:nt)=0;ver2tr(1:nt)=0;exb(1:nx,1)=0;exe(1:nx,1)=0; nf=10;df=1e8;jj=sqrt(-1); gver1in(1:nf)=0;gver1re(1:nf)=0;gver2tr(1:nf)=0; s11(1:nf)=0;s21(1:nf)=0;for n=1:nt for i=1:nx %二維FDTD過程 for k=1:nz hy(i,k)=hy(i,k)+cm*(ez(i+

30、1,k)-ez(i,k)*dz/dx+ex(i,k)-ex(i,k+1); end end if n<=300 ex(1:12,1)=exp(-(n*dt-t0)2/t2); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; %設置擋板 for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy(i,k)-hy(i-1,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; else %高斯波完全進入，入射端改為吸收邊界條件 ex