1、二維TE波FDTD法PML吸收邊界 摘要：由于計算空間有限，一段時間后波面會被反射回來，這是我們不愿看到的。為了避免反射波的存在，需要在計算邊界設置適當的吸收邊界。目前最為有效的吸收邊界是PML層。關鍵詞：二維； TE波；PML層1引言時域有限差分(F D T D )方法已逐漸成為解決電磁散射和電磁波傳輸等問題的有力工具。研究F D T D 方法的核心問題是尋求一種理想的吸收邊界, 使截斷面反射最小。1 9 9 4 年J.P.Berenger 提出了“ 完全匹配層(P M L ) ”這種新邊界完全匹配層技術在計算電磁學領域特別是在FEM中得到廣泛的應用是在Sacks995年提出的單軸各向異性介

2、質PLM以后的事，此后的發展證明，作為一種新型的邊界截斷策略，不管是與FDTD的結合，還是用在FEM中，PLM技術都獲得了巨大成功。2計算模型有一個長和寬都為ke=80的矩形介質，在它的周圍圍上厚度為kp=8的PLM層吸收邊界（如圖一所示），請編寫程序驗證電磁波在經過PLM層時完全被吸收(即電磁波沒有反射)。圖一模型圖 圖二左為程序運行結果Ey圖右為Hzy圖圖三為程序運行結果hz圖3結論：本文的算例證明了PLM作為一種邊界處理的有效性，PLM技術也在有關鍵參數的選取優化問題，這需要研究大量的算例來推出其規律，因而還有很多的工作有待研究。參考文獻：1z S Sacks，D M Kingsland

3、，R Lee and Jin-Fa LeeAperfectly matched anisotropic absorber for use as anabsorbing boundary conditionJIEEE TransAntennasPropagat，1995，43(12)：146014632A Mitchell，T Aberle，D M Kokotoff and M WAustimAn anisotropic PML for use with biaxial mediaDIEEE TransMicrowave Theory Tech，1999，47(3)：3743773班永靈填充各向

4、異性介質的背腔式微帶天線的矢量有限元法分析D北大學碩士論文，2003年5月4班永靈(1978一)，男，河南人，電子科技大學博士研究生，主要研究方向有限元法(特剮是高階矢量有限元)、快速掃頻技術以及并行計算。5周樂柱(1944一)，男，貴州人，北京大學教授、博士生導師，CIE會士，IEEE高級會員，主要研究方向：計算電磁學及其應用(散射、天線、微波器件)，通信中的電磁場問題等。6聶在平(1946一)，男，陜西人，電子科技大學教授，博士生導師，現任電子科技大學副校長，主要研究方向：計算電磁學、電磁散射與逆散射、非均勻介質中的場與渡、移動通信中智能天線技術等。附錄：二維TE波FDTD法PML吸收邊界

5、程序clc,clearke=80;kp=8; ex=ones(ke+2*kp+1,ke+2*kp+1);e=ones(ke+2*kp+1,ke+2*kp+1);ey=ones(ke+2*kp+1,ke+2*kp+1);hz=ones(ke+2*kp+1,ke+2*kp+1);hzx=ones(ke+2*kp+1,ke+2*kp+1);hzy=ones(ke+2*kp+1,ke+2*kp+1); kc=(ke+2*kp)/2;e0=8.85*10-12;mu0=4*3.14*10-7;pi=3.14;c0=3e8;f=1e9;lambda=c0/f;dx=lambda/10;dy=lambda/

6、10;dt=(dx/2)/c0;sinxm=0.333;sinmxm=sigxm*(mu0/e0);nn=3; T=0;for n=1:1000; T=T+1; hz=hzx+hzy; e=ex+ey; hz(kc,kc)=4*sin(2*pi*f*T*dt); for i=1:(2*kp+ke+1); for j=(kp+1):(ke+kp+1); ex(i,j)=ex(i,j)+dt/(e0*dx)*(hz(i,j)-hz(i,j-1); end end for i=1:(ke+2*kp+1); for j=2:kp; sinx=sinxm*(kp+2-j)/kp)nn*2*e0/dt;

7、ex(i,j)=ex(i,j)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0)/(dy*sigx)*sqrt(e0/mu0)*(hz(i,j)-hz(i,j-1); end end for i=1:(ke+2*kp+1); for j=(ke+kp+2):(ke+2*kp+1); sinx=sinxm*(j-ke-kp)/kp)nn*2*e0/dt; ex(i,j)=ex(i,j)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0)/(dy*sigx)*sqrt(e0/mu0)*(hz(i,j)-hz(i,j-1); end end for i=

8、(kp+1):(ke+kp+1); for j=1:(2*kp+ke+1); ey(i,j)=ey(i,j)-dt/(e0*dx)*(hz(i,j)-hz(i-1,j);endend for i=2:kp; for j=1:(ke+2*kp+1); sinx=sinxm*(kp+1-i)/kp)nn*2*e0/dt; ey(i,j)=ey(i,j)*exp(-sigx*dt/e0)-(1-exp(-sigx*dt/e0)/(dx*sigx)*sqrt(e0/mu0)*(hz(i,j)-hz(i-1,j); endend for i=(ke+kp+2):(ke+2*kp+1); for j=1:

9、(ke+2*kp+1); sinx=sinxm*(i-ke-kp-1)/kp)nn*2*e0/dt; ey(i,j)=ey(i,j)*exp(-sinx*dt/e0)-(1-exp(-sigx*dt/e0)/(dx*sigx)*sqrt(e0/mu0)*(hz(i,j)-hz(i-1,j); endend%磁場循環開始 for i=1:(ke+2*kp+1); for j=1:kp; sin(mx)=sinmxm*(kp-j+1)/kp)nn*2*mu0/dt; hzy(i,j)=hzy(i,j)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0)/(dx*sig

10、mx)*sqrt(mu0/e0)*(ex(i,j+1)-ex(i,j); end end for i=1:(2*kp+ke+1); for j=(kp+1):(kp+ke); hzy(i,j)=hzy(i,j)+dt/(mu0*dy)*(ex(i,j+1)-ex(i,j); endend for i=1:(ke+2*kp+1); for j=(kp+ke+1):(ke+2*kp); sigmx=sigmxm*(j-kp-ke)/kp)nn*2*mu0/dt; hzy(i,j)=hzy(i,j)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0)/(dx*sigmx

11、)*sqrt(mu0/e0)*(ex(i,j+1)-ex(i,j); endend for i=1:kp; for j=1:(ke+2*kp+1); sigmx=sigmxm*(kp-i+1)/kp)nn*2*mu0/dt; hzx(i,j)=hzx(i,j)*exp(-sigmx*dt/mu0)-(1-exp(-sigmx*dt/mu0)/(dx*sigmx)*sqrt(mu0/e0)*(ey(i+1,j)-ey(i,j); endend for i=(kp+1):(kp+ke); for j=1:(2*kp+ke+1); hzx(i,j)=hzx(i,j)+dt/(mu0*dy)*(ey(

12、i+1,j)-ey(i,j); endend for i=(ke+kp+1):(ke+2*kp); for j=1:(ke+2*kp+1); sinmx=sinmxm*(i-kp-ke)/kp)nn*2*mu0/dt; hzx(i,j)=hzx(i,j)*exp(-sigmx*dt/mu0)-(1-exp(-sigmx*dt/mu0)/(dx*sigmx)*sqrt(mu0/e0)*(ey(i+1,j)-ey(i,j); endend ey(1,:)=0;ex(:,1)=0;hzx(2*kp+ke+1,:)=0;hzy(:,2*kp+ke+1)=0; plot(hz(:,kc)axis(46,ke+2*kp,45,50) figure(1)subplot(1,2,1)mesh(e