FDTD二維圓柱散射.ppt

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

基于fdtd的PML的TM波散射和mur邊界的TE波散射,組員：樊家偉、黃登祥、袁一粟、 江亞男、陳永煒,我們的任務(wù),1、沒有吸收邊界條件下的TM波同軸圓柱散射； 2、PML吸收邊界條件下的TM波同軸圓柱散射； 3、Mur吸收邊界條件下的TE波圓柱散射。,重寫Maxswell方程組,二維Maxswell方程組標量方程,1、二維FDTD基本原理,對時間和空間差分后，迭代公式,完美匹配層（Perfectly matched Layer, PML）是由Berenger 提出的，使用最為靈活、廣泛的一種ABC（Absorbing Boundary Conditions）。其基本原理是：電磁波的反射量由兩個介質(zhì)的波阻抗決定：,其中,2、完全匹配層基本原理,,2.1 在X 方向上實現(xiàn)PML （僅僅保留與X 方向有關(guān)的ε F* 、 μ F* ）,將ε F* 、 μ F*帶入到左式,注意到： H x方向上的磁導(dǎo)率與 H y方向上的磁導(dǎo)率互成倒數(shù)。因此，滿足 了PML 的第二個條件。,(1),(2),(3),,時域,對(1)式左邊進行差分,,其中,,其中,同理,(2)式,,其中,(3)式,引進輔助參數(shù),隨著電磁波進入到 PML，該參數(shù)是增大的。,i=1,2,…, length_pml,注意到參變量 i/length_pml的變化范圍是從0 到1，而權(quán)值0.333 是保持穩(wěn)定狀態(tài)的最大值,,,,2.2 在y 方向上實現(xiàn)PML,(1),(2),(3),(1),,其中,(2),,其中,(3),源程序 clear all；clc; %設(shè)置網(wǎng)格數(shù) IE = 101; % x 方向網(wǎng)格100,實際有101個點 JE = 101; % y 方向網(wǎng)格多少 %設(shè)置圓柱中心 ic = round(IE / 2); jc = round(JE / 2); %總場區(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; %時域網(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); %電場通量D hx = dz; %x 方向磁場強度 hy = dz; %y 方向磁場強度 ihy = dz; %用于中間變量ihy，是用來計算磁場強度的，是curl_e 的積分 ihx = dz; %用于中間變量ihy，是用來計算磁場強度的，是curl_e 的積分 ga = ones(IE,JE); %電場強度與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); %由于做了歸一化，同時在真空中，因此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 個數(shù)，即PML 有多少個單元網(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ù)可以得到空間和時間上的傳輸情況,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); %計算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;,% 激勵元 在此采用正弦波，頻率為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,%計算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,%計算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; %修正場區(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,%計算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; %修正場區(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,%畫圖 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吸收邊界條件,考慮齊次波動方程，在二維情形為①： 其平面波解為： 其中 設(shè)x=0平面為截斷面，如圖所示，x=0右側(cè)區(qū) 域同時存在入射波和反射波疊加，此區(qū)域中有： ② 式中設(shè) 左行波 右行波,Engquist-Majda吸收邊界條件,②帶入①得 定義微分算子 注：保留對x的導(dǎo)數(shù) 算子L做因式分解 左行波算子 右行波算子,,,Engquist-Majda吸收邊界條件,將左行波算子作用在平面波上得 若在 截斷邊界處設(shè)置條件 ，相當于使截斷截面處右行波（反射波）成分等于零。 將算子具體表示代入上式得 作頻域到時域的轉(zhuǎn)換,,Engquist-Majda吸收邊界條件,Engquist-Majda吸收邊界條件 （截斷邊界位于區(qū)域左側(cè)） （截斷邊界位于區(qū)域右側(cè)） 通過波動方程的因子分解獲得直角坐標系下FDTD吸收邊界的單向波動方程,,,一階和二階近似吸收邊界,算子中含有根式運算，限制了數(shù)值實現(xiàn)。做近似處理： 利用Taylor級數(shù)展開 重寫左行波算子 保留級數(shù)第一項做近似 這就是x=0平面作為區(qū)域左側(cè)界面不產(chǎn)生反射波的一階近似吸收邊界條件。,,一階和二階近似吸收邊界,保留級數(shù)至第二項 這就是x=0平面作為區(qū)域左側(cè)界面不產(chǎn)生反射波的二階近似吸收邊界條件。 二階近似比一階近似吸收邊界條件有所改善，殘留反射波小。 Mur對表中的吸收邊界條件引入了一種簡單有效的差分數(shù)值算法，即對時間和空間的偏微分取二階中心差分近似，將單向波方程離散化，便形成了Mur的一階二階吸收邊界條件，其總體虛假反射在1%~5%，能滿足許多工程需求。,,Mur 吸收邊界條件,對于二維電磁場問題，Mur指出二階近似吸收邊界可降低為只含E,H分量的一階導(dǎo)數(shù)，從而使數(shù)值計算簡化。 TM波 二維直角坐標TM波,,,將上式對t積分，并設(shè)初始場為0,Mur 吸收邊界條件,對于TE波，同理 Mur二階近似吸收邊界條件比一階近似多出一項一階導(dǎo)數(shù)。 Mur吸收邊界條件的FDTD離散式 先看TM波的一階近似 ③,Mur 吸收邊界條件,將上式在 作離散 ④ 利用下式線性插值關(guān)系消去,⑤,Mur 吸收邊界條件,將上式⑤帶入④，再帶入③得 整理后得到TM波邊界條件的一階近似 對照Yee元胞圖可知位于左截斷邊界上的 節(jié)點值是用區(qū)域內(nèi)部節(jié)點值及前一時刻邊界上的節(jié)點值來表示，不涉及截斷邊界以為的場量。,Mur 吸收邊界條件,對于TM波的Mur吸收邊界條件二階形式，比一階多一項 離散式 線性插值 整理后得吸收邊界的離散式：,Mur 吸收邊界條件,TE波一階二階吸收邊界條件可類似推導(dǎo)得到,Mur 吸收邊界條件,對于二維TM波情況，電磁場除了 還有 、 。 由圖可知，在用FDTD計算邊界處的TM波元胞 時并不會涉及截斷邊界以外E或H的節(jié)點。 只有 涉及截斷邊界外側(cè)的H節(jié)點。因此只需給出邊界處切向場分量 的吸收邊界條件。 對TE波同理只需 的 邊界。,Mur 吸收邊界條件,以上只討論了x=0的左側(cè)界面的吸收邊界條件。對于其余三邊有相似結(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)的幅度隨著相對介電常數(shù)不同而呈現(xiàn)不同程度的衰減,TM波在圓柱體內(nèi)的幅度隨著電導(dǎo)率的不同而呈現(xiàn)不同程度的衰減,,沒有吸收邊界時的散射情況：,,PML吸收邊界條件下的散射情況：,