時域有限差分方法林志立實用教案

1、電磁學基本(jbn)方程 麥克斯韋(mi k s wi)方程組( Maxwell Equations)James Clerk Maxwell (18311879)(法拉第感應(gnyng)定律)(安培環(huán)路定律)(高斯定律電場)(高斯定律磁場)l物質(zhì)本構(gòu)關(guān)系 (Constitutive Relations)（歐姆定律）（電極化）（磁化）JE 一切電磁場和電磁波問題均可由以上方程，以及各類具體的邊界條件所決定！第1頁/共28頁第一頁，共28頁。麥克斯韋(mi k s wi)方程組中的運算符散度（Divergence)旋度（Curl)連續(xù)函數(shù)的偏微分運算(yn sun)第2頁/共28頁第二頁，共28

2、頁。設有一連續(xù)函數(shù) , 現(xiàn)欲求(y qi)。FDTD的基本(jbn)思想( )f x( )f x( )fx連續(xù)偏微分的有限(yuxin)階近似時域和空間域的離散化()()( )2 f xxf xxfxx二階中心差分近似表達式：當越小時，上式的近似程度越高。x2()()( )( ).26 f xxf xxxfxfxx實際上：第3頁/共28頁第三頁，共28頁。FDTD空間(kngjin)域的離散化(1)空間(kngjin)域的分割離散化ExEx分量的空分量的空間間(kngjin)(kngjin)離散分布離散分布圖圖YeeYee元胞（元胞（ x, x, y, y, z)z)H Hx x分量的空分量的

3、空間間離散分布離散分布圖圖節(jié)節(jié)點點第4頁/共28頁第四頁，共28頁。FDTD空間(kngjin)域的離散化1()yzxzEHEtyx YEE YEE 元胞元胞各各電電磁磁場場分量分量(fn ling)(fn ling)在元胞中的位在元胞中的位置置K.S. Yee, “Numerical solution of initial boundary value problems involving Maxwells equations in isotropic media,”IEEE Trans. Antennas Propagat., vol. 14, 1966, pp. 302-307.例如例如

4、(lr)(lr)：H Hz z為為ExEx和和EyEy所所環(huán)繞環(huán)繞。第5頁/共28頁第五頁，共28頁。FDTD空間(kngjin)偏微分的近似1()yzxzEHEtyx ( ,1)( , )2()2xxxEEi jEi jyy(1, )( , )2()2yyyEEijEi jxx以HzHz為例：類似地，可實現(xiàn)各電磁場分量(fn ling)的空間偏微分計算。第6頁/共28頁第六頁，共28頁。FDTD時間(shjin)偏微分的近似1()yzxzEHEtyx 以HzHz為例：t=(n+1/2)t=(n+1/2) t tt=nt=n t tt=(n-1/2)t=(n-1/2) t t時間時間(shji

5、n)(shjin)上的推移上的推移(1)(1)電場電場(din ch(din chng)ng)在在時間時間上上取整取整數(shù)數(shù)倍的倍的t; t; t=n t=n * *t;t;(2)(2)磁磁場場在在時間時間上取上取( (整整數(shù)數(shù)1/2)1/2)倍的倍的t;t; t=(n +1/2) t=(n +1/2)* *t;t;1/21/21()2()2nnnnyzzxzEHHEtyx 第7頁/共28頁第七頁，共28頁。麥克斯韋方程(fngchng)的離散化近似1()yzxzEHEtyx 以HzHz為例：1/21/2( , )( , )(1, )( , )1( ,1)( , )() nnnnnnzzyyxx

6、zHi jHi jtEijEi jEi jEi jyx 1/21/2( , )( , )(1, )( , )( ,1)( , )() nnnnnnzzyyxxzHi jHi jEijEi jt Ei jEi jyx 上式即為Hz的更新方程，由前一時刻的磁場和前半時刻的臨近空間(kngjin)格點的電場即可求出最新時刻的磁場。第8頁/共28頁第八頁，共28頁。麥克斯韋方程(fngchng)的離散化近似采取類似的步驟(bzhu)(bzhu)，可以推導出其它場量的更新表達式：11/21/21/21/2( , )( , )( , )(1, )( , )( ,1)() zzynnnnyxxznnEi j

7、Ei jHi jHijtHi jHi jxy 例如(lr)(lr)，對于EzEz：1()yzxzHEHtxy 第9頁/共28頁第九頁，共28頁。FDTD的離散參數(shù)(cnsh)的選擇元胞尺寸元胞尺寸(ch cun)：邊長小于最短波長的：邊長小于最短波長的1/10, 以減小數(shù)值色散。以減小數(shù)值色散。maxmaxmaxmin1, and 10 xyz 222022, ,sin ()sin ()22()()22x y zktct數(shù)值(shz)色散方程：222220()xyzckkk理想色散方程：要求：要求：02k02tmax1 10 tf例如，取例如，取例如，取例如，取第10頁/共28頁第十頁，共28

8、頁。FDTD的離散(lsn)參數(shù)的穩(wěn)定性條件 時間時間(shjin)(shjin)步步長長：Courant Courant 穩(wěn)穩(wěn)定性定性條條件件2221111 tcxyz對對于非色散于非色散(ssn)(ssn)介介質(zhì)質(zhì)，時間時間步步長長不能大于以下表不能大于以下表達達式：式：21/21/222, ,04sin (/2)()( )0,rx y zkZZZct（von Neumann method ）為為了保持了保持穩(wěn)穩(wěn)定性，定性，該該方程的所有解的模必方程的所有解的模必須須小于小于1 1。域數(shù)值色散方程：第11頁/共28頁第十一頁，共28頁。介質(zhì)電磁(dinc)參量的設定( , )( , )1(

9、)yzxzHEHxi ji jty ( , )( , )1()yzxzi jiyjEHEtx 不同的元胞的電磁(dinc)參量應設置為所在空間所代表的介質(zhì)的介電常數(shù)和磁導率。長方體，物質(zhì)I長方體，物質(zhì)II球體，物質(zhì)III空氣場量與介質(zhì)參數(shù)場量與介質(zhì)參數(shù)(cnsh)要對應要對應第12頁/共28頁第十二頁，共28頁。色散(ssn)介質(zhì)的FDTD模擬202200()( )2srj 以以Lorentz 介質(zhì)介質(zhì)(jizh)為例：為例：122012012122012012ee() ,eejjrjjaa Za Zaaabb Zb ZbbbDE1201210012()nnnnnnbbbaaaDDDEEE更新

10、更新(gngxn)方程：方程：近似，近似，求系數(shù)求系數(shù)Z-變換0()()()rDE0( )( )* ( )rtttDE02220002200()( )2( )sin() ( )tsrttet U t 時時域：域：頻頻域：域：第13頁/共28頁第十三頁，共28頁。色散介質(zhì)(jizh)的FDTD模擬模擬模擬Lorentz色散色散(ssn)介質(zhì)的不同方法：介質(zhì)的不同方法：012012,a a a b b bMSE approach :( )( )rrError 12010012120()/ ()/nnnnnnbbbaaaEDDDEE第14頁/共28頁第十四頁，共28頁。FDTD編程流程(lichng

11、)主循主循環(huán)環(huán)(xnhun)(xnhun)初始化初始化輸輸出出(shch)(shch)結(jié)結(jié)果果第15頁/共28頁第十五頁，共28頁。編程舉例(j l)1：一維FDTD問題x x基本基本(jbn)(jbn)旋度方程：旋度方程：X向電導率X向磁電導率X向電流(dinli)X向磁流第16頁/共28頁第十六頁，共28頁。編程舉例編程舉例(j l)1：一維：一維FDTD問題問題Matlab程序代碼：% Define initial constantseps_0 = 8.854187817e-12; % permittivity of free space mu_0 = 4*pi*1e-7; % perm

12、eability of free space c = 1/sqrt(mu_0*eps_0); % speed of light % Define problem geometry and parametersdomain_size = 1; % 1D problem space length in metersdx = 1e-3; % cell size in meters dt = 3e-12; % duration of time step in seconds number_of_time_steps = 2000; % number of iterations nx = round(d

13、omain_size/dx); % number of cells in 1D problem spacesource_position = 0.5; % position of the current source Jz 1.定義(dngy)物理常量2.定義問題的參量和結(jié)構(gòu)(jigu)尺寸第17頁/共28頁第十七頁，共28頁。% Initialize field and material arraysCeze = zeros(nx+1,1);Cezhy = zeros(nx+1,1);Cezj = zeros(nx+1,1);Ez = zeros(nx+1,1);Jz = zeros(nx+

14、1,1);eps_r_z = ones (nx+1,1); % free spacesigma_e_z = zeros(nx+1,1); % free space Chyh = zeros(nx,1);Chyez = zeros(nx,1);Chym = zeros(nx,1);Hy = zeros(nx,1);My = zeros(nx,1);mu_r_y = ones (nx,1); % free spacesigma_m_y = zeros(nx,1); % free space 編程舉例編程舉例(j l)1：一維：一維FDTD問題（續(xù)）問題（續(xù)）3.初始化場量和介質(zhì)(jizh)參量陣列

15、電場電電場電流流(dinli)部分部分磁磁場場磁流部分磁流部分第18頁/共28頁第十八頁，共28頁。% Calculate FDTD updating coefficients Ceze = (2 * eps_r_z * eps_0 - dt * sigma_e_z) . ./(2 * eps_r_z * eps_0 + dt * sigma_e_z); Cezhy = (2 * dt / dx) . ./(2 * eps_r_z * eps_0 + dt * sigma_e_z); Cezj = (-2 * dt) . ./(2 * eps_r_z * eps_0 + dt * sigma_

16、e_z); Chyh = (2 * mu_r_y * mu_0 - dt * sigma_m_y) . ./(2 * mu_r_y * mu_0 + dt * sigma_m_y); Chyez = (2 * dt / dx) . ./(2 * mu_r_y * mu_0 + dt * sigma_m_y); Chym = (-2 * dt) . ./(2 * mu_r_y * mu_0 + dt * sigma_m_y); 編程舉例編程舉例(j l)1：一維：一維FDTD問題（續(xù)）問題（續(xù)）4.計算更新方程(fngchng)系數(shù)電場電場(din chng)部分部分磁場部分第19頁/共28頁第

17、十九頁，共28頁。% Define the Gaussian source waveform time = dt*0:number_of_time_steps-1.;Jz_waveform = exp(-(time-2e-10)/5e-11).2);source_position_index = round(nx*source_position/domain_size)+1; % Subroutine to initialize plotting initialize_plotting_parameters; 編程舉例編程舉例(j l)1：一維：一維FDTD問題（續(xù)）問題（續(xù)）5.定義(dng

18、y)場源源波形(b xn)為高斯型6.作圖初始化Ez_positions = 0:nx*dx;Hy_positions = (0:nx-1+0.5)*dx;v = 0 -0.1 -0.1; 0 -0.1 0.1; 0 0.1 0.1; 0 0.1 -0.1; . 1 -0.1 -0.1; 1 -0.1 0.1; 1 0.1 0.1; 1 0.1 -0.1;f = 1 2 3 4; 5 6 7 8;axis(0 1 -0.2 0.2 -0.2 0.2);lez = line(Ez_positions,Ez*0,Ez,Color,b,LineWidth,1.5);lhy = line(Hy_pos

19、itions,377*Hy,Hy*0,Color,r, . LineWidth,1.5,linestyle,-.);set(gca,fontsize,12,FontWeight,bold);axis square;legend(E_z, H_y times 377,Location,NorthEast);xlabel(x m);ylabel(A/m);zlabel(V/m);grid on;p = patch(vertices, v, faces, f, facecolor, g, facealpha,0.2);text(0,1,1.1,PEC,horizontalalignment,cent

20、er,fontweight,bold);text(1,1,1.1,PEC,horizontalalignment,center,fontweight,bold);第20頁/共28頁第二十頁，共28頁。% FDTD loopfor time_step = 1:number_of_time_steps % Update Jz for the current time step Jz(source_position_index) = Jz_waveform(time_step); % Update magnetic field Hy(1:nx) = Chyh(1:nx) .* Hy(1:nx) .

21、+ Chyez(1:nx) .* (Ez(2:nx+1) - Ez(1:nx) . + Chym(1:nx) .* My(1:nx); % Update electric field Ez(2:nx) = Ceze (2:nx) .* Ez(2:nx) . + Cezhy(2:nx) .* (Hy(2:nx) - Hy(1:nx-1) . + Cezj(2:nx) .* Jz(2:nx); Ez(1) = 0; % Apply PEC boundary condition at x = 0 m Ez(nx+1) = 0; % Apply PEC boundary condition at x

22、= 1 m % Subroutine to plot the current state of the fields plot_fields;end編程舉例(j l)1：一維FDTD問題（續(xù)）7.主循環(huán)(xnhun)程序引入電流(dinli)源更新磁場更新電場設置完美電導層跟新圖面第21頁/共28頁第二十一頁，共28頁。% subroutine used to plot 1D transient fields delete(lez);delete(lhy);lez = line(Ez_positions,Ez*0,Ez,Color,b,LineWidth,1.5);lhy = line(Hy_

23、positions,377*Hy,Hy*0,Color,r, . LineWidth,1.5,linestyle,-.);ts = num2str(time_step);ti = num2str(dt*time_step*1e9);title(time step = ts , time = ti ns);drawnow;編程舉例(j l)1：一維FDTD問題（續(xù)）7.更新(gngxn)畫面運運行行(ynxng)(ynxng)程序程序: :fdtd_1d_code.mfdtd_1d_code.m第22頁/共28頁第二十二頁，共28頁。開放(kifng)模擬空間的模擬完美匹配完美匹配(ppi)層（層（Perfectly Matched Layer)的概念的