FDTD二維圓柱散射.ppt

《FDTD二維圓柱散射.ppt》由會(huì)員分享，可在線閱讀，更多相關(guān)《FDTD二維圓柱散射.ppt（53頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

基于fdtd的PML的TM波散射和mur邊界的TE波散射,組員：樊家偉、黃登祥、袁一粟、 江亞男、陳永煒,我們的任務(wù),1、沒(méi)有吸收邊界條件下的TM波同軸圓柱散射； 2、PML吸收邊界條件下的TM波同軸圓柱散射； 3、Mur吸收邊界條件下的TE波圓柱散射。,重寫(xiě)Maxswell方程組,二維Maxswell方程組標(biāo)量方程,1、二維FDTD基本原理,對(duì)時(shí)間和空間差分后，迭代公式,完美匹配層（Perfectly matched Layer, PML）是由Berenger 提出的，使用最為靈活、廣泛的一種ABC（Absorbing Boundary Conditions）。其基本原理是：電磁波的反射量由兩個(gè)介質(zhì)的波阻抗決定：,其中,2、完全匹配層基本原理,,2.1 在X 方向上實(shí)現(xiàn)PML （僅僅保留與X 方向有關(guān)的ε F* 、 μ F* ）,將ε F* 、 μ F*帶入到左式,注意到： H x方向上的磁導(dǎo)率與 H y方向上的磁導(dǎo)率互成倒數(shù)。因此，滿足 了PML 的第二個(gè)條件。,(1),(2),(3),,時(shí)域,對(duì)(1)式左邊進(jìn)行差分,,其中,,其中,同理,(2)式,,其中,(3)式,引進(jìn)輔助參數(shù),隨著電磁波進(jìn)入到 PML，該參數(shù)是增大的。,i=1,2,…, length_pml,注意到參變量 i/length_pml的變化范圍是從0 到1，而權(quán)值0.333 是保持穩(wěn)定狀態(tài)的最大值,,,,2.2 在y 方向上實(shí)現(xiàn)PML,(1),(2),(3),(1),,其中,(2),,其中,(3),源程序 clear all；clc; %設(shè)置網(wǎng)格數(shù) IE = 101; % x 方向網(wǎng)格100,實(shí)際有101個(gè)點(diǎn) JE = 101; % y 方向網(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ù)PI c=3e8;,%初始化系統(tǒng)變量 dz = zeros(IE,JE); %電場(chǎng)通量D hx = dz; %x 方向磁場(chǎng)強(qiáng)度 hy = dz; %y 方向磁場(chǎng)強(qiáng)度 ihy = dz; %用于中間變量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;ez_inc_high_m1=0; ez_inc_high_m2=0; radius = 15; epsilon=10;%越大衰減越大 sigma=0.3;%越大衰減越大 for j=ja:jb for i=ia:ib xdist=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; end end end,%內(nèi)圓 radius1 = 10; epsilon1=1000; sigma1=10; for j=ja:jb for i=ia:ib xdist=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; end end end,%輸入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 = 1:npml xnum = npml -i + 1; %從npml 到0 xd = npml; xxn = xnum / xd; %輔助變量xxn，從1 到0 xn = 0.33 * xxn^3; %成立方衰減 gi2(i) = 1 / ( 1 + xn ); gi2(IE-i+1) = 1 / ( 1 + xn ); gi3(i) = (1 - xn) / (1 + xn); gi3(IE-i+1) = ( 1 - xn ) / ( 1 + xn ); xxn = ( xnum - 0.5 ) / xd; if(xxn0) break; end,xn = 0.33 * xxn^3; 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+1 xnum = npml - j + 1; xd = npml; xxn = xnum / xd; xn = 0.33 * xxn^3; gj2(j) = 1 / ( 1 + xn ); gj2(JE+1-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) break end xn = 0.33 * xxn^3; 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 + xn ); 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(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); %計(jì)算D for j = 2:JE 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;,% 激勵(lì)元 在此采用正弦波，頻率為1.5G，可選為高斯脈沖 % hard source %pulse_choice = input('Enter 1 for Gaussian, 2 for sine wave pulse: '); %if pulse_choice == 1 pulse = 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ì)算E for j = 1:JE for i = 1:IE ez(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 PML for j=1:JE ez(1,j) = 0; ez(IE,j) = 0; end; for i = 1:IE ez(i,1) = 0; ez(i,JE) = 0; end; for j=1:JE-1 hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1)); end for j=1:JE for i=1:IE real_pt(i,j)=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); end end,%計(jì)算Hx for j = 1:JE-1 for i = 1:IE curl_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:ib hx(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ì)算Hy for j = 1:JE for i = 1:IE -1 curl_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:jb hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j); hy(ib,j)=hy(ib,j)+0.5*ez_inc(j); end,%畫(huà)圖 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_in^2+imag_in^2); for j=ja:jb if(ga(ic,j)1) amp(j)=sqrt(real_pt(ic,j)^2+imag_pt(ic,j)^2)/amp_in; end end figure(2); plot(ja+26:jb-30,amp(ja+26:jb-30)); title('電磁波的幅度變化');,3、Mur 吸收邊界原理,在此之前先了解一下 Engquist-Majda吸收邊界條件,Engquist-Majda吸收邊界條件,考慮齊次波動(dòng)方程，在二維情形為①： 其平面波解為： 其中 設(shè)x=0平面為截?cái)嗝妫鐖D所示，x=0右側(cè)區(qū) 域同時(shí)存在入射波和反射波疊加，此區(qū)域中有： ② 式中設(shè) 左行波 右行波,Engquist-Majda吸收邊界條件,②帶入①得 定義微分算子 注：保留對(duì)x的導(dǎo)數(shù) 算子L做因式分解 左行波算子 右行波算子,,,Engquist-Majda吸收邊界條件,將左行波算子作用在平面波上得 若在 截?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ù)展開(kāi) 重寫(xiě)左行波算子 保留級(jí)數(shù)第一項(xiàng)做近似 這就是x=0平面作為區(qū)域左側(cè)界面不產(chǎn)生反射波的一階近似吸收邊界條件。,,一階和二階近似吸收邊界,保留級(jí)數(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)為0,Mur 吸收邊界條件,對(duì)于TE波，同理 Mur二階近似吸收邊界條件比一階近似多出一項(xiàng)一階導(dǎo)數(shù)。 Mur吸收邊界條件的FDTD離散式 先看TM波的一階近似 ③,Mur 吸收邊界條件,將上式在 作離散 ④ 利用下式線性插值關(guān)系消去,⑤,Mur 吸收邊界條件,將上式⑤帶入④，再帶入③得 整理后得到TM波邊界條件的一階近似 對(duì)照Yee元胞圖可知位于左截?cái)噙吔缟系? 節(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波一階二階吸收邊界條件可類(lèi)似推導(dǎo)得到,Mur 吸收邊界條件,對(duì)于二維TM波情況，電磁場(chǎng)除了 還有 、 。 由圖可知，在用FDTD計(jì)算邊界處的TM波元胞 時(shí)并不會(huì)涉及截?cái)噙吔缫酝釫或H的節(jié)點(diǎn)。 只有 涉及截?cái)噙吔缤鈧?cè)的H節(jié)點(diǎn)。因此只需給出邊界處切向場(chǎng)分量 的吸收邊界條件。 對(duì)TE波同理只需 的 邊界。,Mur 吸收邊界條件,以上只討論了x=0的左側(cè)界面的吸收邊界條件。對(duì)于其余三邊有相似結(jié)果。,Mur程序 網(wǎng)格劃分,% Define the FDTD grid s_smaller = 1; t_smaller = sqrt(2); del_s = lambda/20/s_smaller del_t = del_s/c/sqrt(2)/t_smaller dim_s = 5; s_range = dim_s * lambda; s_range = ceil(s_range / del_s) * del_s; s_cells = s_range / del_s; cells = s_cells; nodes = cells+1;,主循環(huán),while (done~=0), % Store previous values of E Ex_o = Ex; Ey_o = Ey; Hz_o = Hz;,Mur邊界,% Update boundary conditions (Ex , Ey) Ex(1,1:nodes) = Ex_o(2,1:nodes) + c_MUR*(Ex(2,1:nodes) - Ex(1,1:nodes)); Ey(1:nodes,nodes) = Ey_o(1:nodes,nodes-1) + c_MUR*(Ey(1:nodes,nodes-1) - Ey(1:nodes,nodes)); Ex(nodes,1:nodes) = Ex_o(nodes-1,1:nodes) + c_MUR*(Ex(nodes-1,1:nodes) - Ex(nodes,1:nodes)); Ey(1:nodes,1) = Ey_o(1:nodes,2) + c_MUR*(Ey(1:nodes,2) - Ey(1:nodes,1));,作圖,% Plot response figure(5); mesh(abs(Ex + i*Ey)); axis([0 nodes 0 nodes 0 1]); end;,TE波MUR效果,5、仿真結(jié)果及分析,TM波在圓柱體內(nèi)的幅度隨著相對(duì)介電常數(shù)不同而呈現(xiàn)不同程度的衰減,TM波在圓柱體內(nèi)的幅度隨著電導(dǎo)率的不同而呈現(xiàn)不同程度的衰減,,沒(méi)有吸收邊界時(shí)的散射情況：,,PML吸收邊界條件下的散射情況：,