FD二維圓柱散射

1、基于fdtd的PML的TM波散射和mur邊界的TE波散射組員：樊家偉、黃登祥、袁一粟、江亞男、陳永煒我們的任務(wù)? 1、沒有吸收邊界條件下的TM波同軸圓柱散射；? 2、PML吸收邊界條件下的TM波同軸圓柱散射；? 3、Mur吸收邊界條件下的TE波圓柱散射。1、二維FDTD基本原理重寫Maxswell方程組?D1? ? H?t?0?0*rD?=? ?E? ? ?H1? ? ? E?t?0?0二維Maxswell方程組標(biāo)量方程H?H?D1?yxz?t?x?y?0?0?D?=?E? ? ? ? ?zz*r?HxE1?z?t?0?0?y?Hy?Ez?t?0?0?x1對(duì)時(shí)間和空間差分后，迭代公式對(duì)時(shí)間和空

2、間差分后，迭代公式1?tnnD?i, j ? Di, j ?Hi ?1/ 2, j ? Hi ?1/ 2, j? ?y?y?x?00n?1/2zn?1/2z1?tnn?Hi, j?1/ 2 ? Hi, j?1/ 2?x?x?y?001?tn?1/2n?1/2Hi, j ? Hi, j ?Ei, j?1 ? Ei, j?z?z?y00n?1xnx?1?tn?1/2n?1/2Hi, j ? Hi, j ?Ei ? 1, ji, j?E?z?z?x00n?1yny?n?1/2n?1/2nnDi, j ? Di, j ? 0.5 Hi?1/ 2, j ? Hi?1/ 2, j? ? ?y?zzynn?

3、 0.5 Hi, j?1/ 2 ? Hi, j?1/ 2?x?xn?1nn?1/2n?1/2Hi, j ? Hi, j ?0.5 Ei,1j? Ei, j? ? ?xxzz?n?1nn?1/2n?1/2Hi, j ? Hi, j ?0.5 Eij?1,? Ei, j? ? ?yyzz?2、完全匹配層基本原理完美匹配層（Perfectly matched Layer, PML）是由Berenger 提出的，使用最為靈活、廣泛的一種ABC（Absorbing Boundary Conditions）。其基本原理是：電磁波的反射量由兩個(gè)介質(zhì)的波阻抗決定：?A?B?A?B其中? ?2.1 在X 方向上

4、實(shí)現(xiàn)PML（僅僅保留與X 方向有關(guān)的F* 、 F* ）?H?H1?y?xj?D?x?x?y?00?*zFz?Ezj?H?xc? ?0?y*xFx將F* 、 F*帶入到左式?H?H() x ?y?xDj?D1? c?z?0?x?y?0?j?(x)?EDzj?H1? ?c?x?0?y0?j?(x)?EDzj?H1? c?y0j?x0?1?Ezj?H?x? c0?x*yFy(1)(2)(3)注意到： H x方向上的磁導(dǎo)率與 H y方向上的磁導(dǎo)率互成倒數(shù)。因此，滿足了PML 的第二個(gè)條件。?H?H(x)?y?xj?D1? c?z?0?x?y?0?j對(duì)(1)式左邊進(jìn)行差分?() x() x ?Dj?D1

5、? j?DD?z?z?z?0?0?j時(shí)域?Dz?D(x)?Dz?t?0n?1/2n?1/2n?1/2n?1/2Di, j ? Di, jDi, j ? Di, j?(i)z? ?z?z? ?z? ?D?t20?D?n?D?(i)?t(i)?t11?1/2? Di, j1? Di, j1? ?z?t?t0?0?2?2n?1/2z?H?H(x)?y?xj?D1? c?z?0?x?y?0?jnn?Di, j ? gii 3Di, j ? gii 20.5 Hi ?1/2, j ? Hi?1/2, j? ? ?yy?nn?Hi, j?1/2 ? Hi, j?1/ 2?xx?n?1/2zn?1/2z其中

6、?t1c0?x21gi2 i ?1?i ?t/?2?D01?i ?t/?2?D0gi3 i ?1?i ?t/?2?D0?(x)?EDzj?H1? c?0y?j?x0?(2)式同理11n?n?1111?n?1n22Hi?, jfi3 i?Hi?, jfii 20.5 Ei?1, jE?i, j?yyz?z?2222?其中?t1c0?x21?1 ?fi2 i?1 ?2 ?1?i?t/?2?D0?2 ?1 ?1?D?i?t/?2?0?2?1 ?fi3?i?2?1?i?1 ?0?t/?2D?2?(x)?EDzj?H1? ?c?x?0?y0?j?1(3)式1?1?1?n?H?i, j? Hurl _ex

7、?i, j?c2?2?2?1 ?n?1/2? fi1?i?IHx?i, j?2?n?1x其中curl _e? Ei, j ? Ei, j?1?n?1/2zn?1/2x11?n?1/2I?ij ,? Iij ,? curl_eHx?22?n?1/2Hx?t?D?i?fi1?i?2?0c0?t? 0.5?x引進(jìn)輔助參數(shù)?t?Dxn?2?0隨著電磁波進(jìn)入到 PML，該參數(shù)是增大的。?ixni.333? 0?ength_ pm l?l3i=1,2, length_pml1gi2?i?1? xn?i?1? xn?i?gi3?i?1? xn?i?1,0.75?1,0.5?0,0.333?fi1?i?注意到

8、參變量 i/length_pml的變化范圍是從0 到1，而權(quán)值0.333 是保持穩(wěn)定狀態(tài)的最大值2.2 在y 方向上實(shí)現(xiàn)PML(1)?H?H(x)() yy?xDDj?1?1?D ? c?z0?x?y?00?j?j? ?() x()y ?EDDzj?H1?1? ?c? ?0 x?j?y0? ?0?j?() x ?() y ?EDDzj?1?1?H ? c?y0?x0?0?j?j?1?1(2)(3)?H?(1)H(x)() yy?xDDj?1?1?D ? c?z0?x?y?00?j?j?n?1/2n?1/2Di, j ? gi3 i gj3 j Di, j? ?z? ?znn? gi2 i gj

9、2 j 0.5 Hi1/2, j ? Hi1/ 2, j?yy?nn?Hi, j?1/2 ? Hi1/ 2?x?, j?x? ?() x() y ?EDDzj?H1?1? ?c(2)? ?0 x?j?y0? ?0?j?111?1 ?n?1?1 ?H?i, j? fj3Hi, j?fj2curl _ex?j?j?22?2?2?2?1?n?1/2? fi1 i Ii, j?Hx?2?n?1x其中curl _e? Ei, j ? Ei, j?1?n?1/2zn?1/2x11?n?1/2I?ij ,? Iij ,? curl_eHx?22?n?1/2Hx?t?D?i?fi1?i?2?0c0?t? 0.

10、5?x?() x ?() y ?EDDzj?1?1?H ? c?y0?x0?0?j?j?1(3)Hn?1y?1?1?n?1?i?, j? fi3?i?Hy?i?, j?2?2?2?1? fi2?i?0.5curl_e?2?1?1 ?n?/12? fi1?i?0.5IHy?i?, j?2?2?其中curl _e? Ei, j ? Ei, j?1?n?1/2xn?1/2x11?n?1/2I?ij , ? Iij , ?curl _eHy?22?n?1/2Hy?源程序?clear all；clc;%設(shè)置網(wǎng)格數(shù)IE = 101; % x 方向網(wǎng)格100,實(shí)際有101個(gè)點(diǎn)JE = 101; % y 方向

11、網(wǎng)格多少%設(shè)置圓柱中心ic = round(IE / 2);jc = round(JE / 2);%總場(chǎng)區(qū)區(qū)域ia=15;ib=IE-ia+1;ja =15;jb=JE-ja+1;%初始化設(shè)置ddx = .01; %空間網(wǎng)格大小，x 方向和y 方向網(wǎng)格大小相同dt = ddx / 6e8; %時(shí)域網(wǎng)格大小epsz = 8.8e-12; %介電常數(shù)epsilon=30; sigma=0.3;pi = 3.13159; %常數(shù)PIc=3e8;?%初始化系統(tǒng)變量dz = zeros(IE,JE); %電場(chǎng)通量Dhx = dz; %x 方向磁場(chǎng)強(qiáng)度hy = dz; %y 方向磁場(chǎng)強(qiáng)度ihy = dz;

12、 %用于中間變量ihy，是用來(lái)計(jì)算磁場(chǎng)強(qiáng)度的，是curl_e 的積分ihx = dz; %用于中間變量ihy，是用來(lái)計(jì)算磁場(chǎng)強(qiáng)度的，是curl_e 的積分ga = ones(IE,JE); %電場(chǎng)強(qiáng)度與D 之間的關(guān)系矩陣，gb=zeros(IE,JE);iz=gb;real_pt=gb;imag_pt=gb;real_in=0;imag_in=0;amp=zeros(JE,1);%由于做了歸一化，同時(shí)在真空中，因此ga 為1?%初始化變量ez_inc=zeros(JE,1);hx_inc=ez_inc;hy_inc=ez_inc;ez_inc_low_m1=0;ez_inc_low_m2=0;

13、ez_inc_high_m1=0;ez_inc_high_m2=0;radius = 15;epsilon=10;%越大衰減越大sigma=0.3;%越大衰減越大for j=ja:jbfor i=ia:ibxdist=ic-i;ydist=jc-j;dist=sqrt(xdist*xdist+ydist*ydist);if(dist=radius)ga(i,j)=1./(epsilon+(sigma*dt/epsz);gb(i,j)=sigma*dt/epsz;endendend?%內(nèi)圓radius1 = 10;epsilon1=1000;sigma1=10;for j=ja:jbfor i=

14、ia:ibxdist=ic-i;ydist=jc-j;dist=sqrt(xdist*xdist+ydist*ydist);if(dist=radius1)ga(i,j)=1./(epsilon1+(sigma1*dt/epsz);gb(i,j)=sigma1*dt/epsz;endendend?%輸入PML Cell 個(gè)數(shù)，即PML 有多少個(gè)單元網(wǎng)格，在此，x 方向和y 方向上的PML 網(wǎng)格相同?npml = input(Please input the number of PML Cell: );?%x 方向用*i*表示，一方向用*j*表示?%x 方向上的PML 參數(shù)設(shè)置?for i =

15、1:npml?xnum = npml -i + 1; %從npml 到0?xd = npml;?xxn = xnum / xd; %輔助變量xxn，從1 到03?xn = 0.33 * xxn3; %成立方衰減?i?gi2(i) = 1 / ( 1 + xn );xni0.333?ength_ pm l?gi2(IE-i+1) = 1 / ( 1 + xn );?l?gi3(i) = (1 - xn) / (1 + xn);?gi3(IE-i+1) = ( 1 - xn ) / ( 1 + xn );?xxn = ( xnum - 0.5 ) / xd;?if(xxn0)?break;?end

16、?xn = 0.33 * xxn3;fil(i) = xn;fil(IE-i+1) = xn;fi2(i) = 1 / ( 1 + xn );fi2(IE-i+1) = 1 / ( 1 + xn );fi3(i) = ( 1-xn ) / ( 1 + xn );fi3(IE-i+1) = ( 1 - xn ) / ( 1 + xn );end;%y 方向上的PML 參數(shù)設(shè)置for j = 1:npml+1xnum = npml - j + 1;xd = npml;xxn = xnum / xd;xn = 0.33 * xxn3;gj2(j) = 1 / ( 1 + xn );gj2(JE+1-

17、j) = 1 / ( 1 + xn );gj3(j) = ( 1 - xn ) / ( 1 + xn );gj3(JE-j+1) = ( 1 - xn ) / ( 1 + xn );xxn = ( xnum - 0.5 ) / xd;?if(xxn0)breakendxn = 0.33 * xxn3;fj1(j) = xn;fj1(JE-j+1) = xn;fj2(j) = 1 / ( 1 + xn );fj2(JE-j+1) = 1 / ( 1 + xn );fj3(j) = ( 1 - xn ) / ( 1 + xn );fj3(JE-j+1) = ( 1 - xn ) / ( 1 + x

18、n );end;%高斯脈沖變量設(shè)置t0 = 40;spread = 10;T = 0;%輸入nsteps 必須為正整數(shù)nsteps = input(Please input the number of nsteps );%設(shè)置循環(huán)次數(shù)，從常數(shù)可以得到空間和時(shí)間上的傳輸情況? for n = 1:nsteps?T = T+1;?for j=2:JE?ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j);?end?real_in=real_in+cos(2*pi*frequency*T*dt)*ez_inc(ja);?imag_in=imag_in-sin(

19、2*pi*frequency*T*dt)*ez_inc(ja);?%邊界條件?ez_inc(1)=ez_inc_low_m2;?ez_inc_low_m2=ez_inc_low_m1;?ez_inc_low_m1=ez_inc(1);?ez_inc(JE)=ez_inc_high_m2;?ez_inc_high_m2=ez_inc_high_m1;?ez_inc_high_m1=ez_inc(JE-1);D?ij ,? gi3 i gj3 j D?ij ,? ?%計(jì)算D?gi2 i gj2 j 0.5 Hi1/2, j ?Hi1/2, j? ?for j = 2:JE?Hij1/2 ?Hij1

20、/2?, ?, ?for i = 2:IE?dz(i,j) = gi3(i) * gj3(j) * dz(i,j) + gi2(i) * gj2(j) * 0.5 * ( hy(i,j) - hy(i-1,j) - hx(i,j) + hx(i,j-1) );?end;?end;n?1/2zn?1/2znynynxnx?% 激勵(lì)元 在此采用正弦波，頻率為1.5G，可選為高斯脈沖?% hard source?%pulse_choice = input(Enter 1 for Gaussian, 2 for sine wave pulse: );?%if pulse_choice = 1?puls

21、e = exp( -0.5 * ( t0-T ) / spread ) 2 ) ;?%clear j;?%pulse=1;%exp(-j*omiga*ddx*T/c);?%dz(IE/2,JE/2) = pulse + dz(IE/2,JE/2);?%else?%pulse = sin( 2*pi*1500*1e6*dt*T);?ez_inc(10) = pulse;?%end;?for i=ia:ib?dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1);?dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);?end?%計(jì)算Efor j = 1:JEfor i

22、= 1:IEez(i,j) = ga(i,j) * (dz(i,j)-iz(i,j);iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j);end;end;%Set the Ez edges to 0, as part of the PMLfor j=1:JEez(1,j) = 0;ez(IE,j) = 0;end;for i = 1:IEez(i,1) = 0;ez(i,JE) = 0;end;for j=1:JE-1hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1);endfor j=1:JEfor i=1:IEreal_pt(i,j)=

23、real_pt(i,j)+cos(2*pi*frequency*T*dt)*ez(i,j);imag_pt(i,j)=imag_pt(i,j)-sin(2*pi*frequency*T*dt)*ez(i,j);endend?%計(jì)算Hxfor j = 1:JE-1for i = 1:IEcurl_e = ez(i,j) - ez(i,j+1);ihx(i,j) = ihx(i,j) + fi1(i) * curl_e;hx(i,j) = fj3(j) * hx(i,j) + fj2(j) * 0.5 * ( curl_e + ihx(i,j) );end; end;%修正場(chǎng)區(qū)域for i=ia:

24、ibhx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja);hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb);end?%計(jì)算Hyfor j = 1:JEfor i = 1:IE -1curl_e = ez(i+1,j) - ez(i,j);ihy(i,j) = ihy(i,j) + fj1(j) * curl_e;hy(i,j) = fi3(i) * hy(i,j) + fi2(i) * 0.5 * ( curl_e + ihy(i,j) );end;end;%修正場(chǎng)區(qū)域for j=ja:jbhy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j

25、);hy(ib,j)=hy(ib,j)+0.5*ez_inc(j);end?%畫圖figure(1);surfc(1:IE,1:JE,ez)%surfc(ia:ib,ja:jb,ez(ia:ib,ja:jb)axis(0 IE 0 JE -0.2 1.)pause(.01);end;amp_in=sqrt(real_in2+imag_in2);for j=ja:jbif(ga(ic,j)=0右側(cè)區(qū)域同時(shí)存在入射波和反射波疊加，此區(qū)域中有：式中設(shè)左行波右行波Engquist-Majda吸收邊界條件?帶入得定義微分算子注：保留對(duì) x的導(dǎo)數(shù)?算子L做因式分解左行波算子右行波算子Engquist-Ma

26、jda吸收邊界條件?將左行波算子作用在平面波上得若在 截?cái)噙吔缣幵O(shè)置條件行波（反射波）成分等于零。，相當(dāng)于使截?cái)嘟孛嫣幱?將算子具體表示代入上式得?作頻域到時(shí)域的轉(zhuǎn)換Engquist-Majda吸收邊界條件?Engquist-Majda吸收邊界條件（截?cái)噙吔缥挥趨^(qū)域左側(cè)）（截?cái)噙吔缥挥趨^(qū)域右側(cè)）?通過(guò)波動(dòng)方程的因子分解獲得直角坐標(biāo)系下FDTD吸收邊界的單向波動(dòng)方程一階和二階近似吸收邊界?算子中含有根式運(yùn)算，限制了數(shù)值實(shí)現(xiàn)。做近似處理：?利用Taylor級(jí)數(shù)展開?重寫左行波算子?保留級(jí)數(shù)第一項(xiàng)做近似這就是x=0平面作為區(qū)域左側(cè)界面不產(chǎn)生反射波的一階近似吸收邊界條件。一階和二階近似吸收邊界?保留級(jí)

27、數(shù)至第二項(xiàng)這就是x=0平面作為區(qū)域左側(cè)界面不產(chǎn)生反射波的二階近似吸收邊界條件。?二階近似比一階近似吸收邊界條件有所改善，殘留反射波小。?Mur對(duì)表中的吸收邊界條件引入了一種簡(jiǎn)單有效的差分?jǐn)?shù)值算法，即對(duì)時(shí)間和空間的偏微分取二階中心差分近似，將單向波方程離散化，便形成了Mur的一階二階吸收邊界條件，其總體虛假反射在1%5%，能滿足許多工程需求。Mur 吸收邊界條件?對(duì)于二維電磁場(chǎng)問(wèn)題，Mur指出二階近似吸收邊界可降低為只含E,H分量的一階導(dǎo)數(shù)，從而使數(shù)值計(jì)算簡(jiǎn)化。?TM波二維直角坐標(biāo)TM波將上式對(duì)t積分，并設(shè)初始場(chǎng)為0Mur 吸收邊界條件?對(duì)于TE波，同理Mur二階近似吸收邊界條件比一階近似多出一

28、項(xiàng)一階導(dǎo)數(shù)。?Mur吸收邊界條件的FDTD離散式?先看TM波的一階近似Mur 吸收邊界條件i?1/2，j)及t? (n?1/ 2)?t作離散? 將上式在(EZ? 利用下式線性插值關(guān)系消去Mur 吸收邊界條件?將上式帶入，再帶入得?整理后得到TM波邊界條件的一階近似?對(duì)照Yee元胞圖可知位于左截?cái)噙吔缟系腅Z節(jié)點(diǎn)值是用區(qū)域內(nèi)部節(jié)點(diǎn)值及前一時(shí)刻邊界上的節(jié)點(diǎn)值來(lái)表示，不涉及截?cái)噙吔缫詾榈膱?chǎng)量。Mur 吸收邊界條件?對(duì)于TM波的Mur吸收邊界條件二階形式，比一階多一項(xiàng)離散式線性插值?整理后得吸收邊界的離散式：Mur 吸收邊界條件?TE波一階二階吸收邊界條件可類似推導(dǎo)得到Mur 吸收邊界條件?對(duì)于二維TM波情況，電磁場(chǎng)除了EZ還有Hx、 。Hy由圖可知，在用FDTD計(jì)算邊界處的TM波元胞時(shí)并不會(huì)涉及截?cái)噙吔缫酝釫或H的節(jié)點(diǎn)。EZ只有涉及截?cái)噙吔缤鈧?cè)的H節(jié)點(diǎn)。因此只需給出邊界處切向場(chǎng)分量的吸收邊界條件。Hz?對(duì)TE波同理只需的邊界。Mur 吸收邊