時域有限差分方法-林志立

1、Beihang University時域有限差分方法林志立北航儀器光電學院光電工程系The Finite-Difference Time-Domain Method for Electromagnetics計算電磁學中的Beihang University電磁學中幾種重要的數(shù)值計算方法l有限差分法 Finite Difference Method靜電場、靜磁場的有限差分法；時域行波的電磁場的時域有限差分法；l有限元法 (Finite Element Method)數(shù)值求解各類獨立的偏微分方程；（電磁學、材料力學、工程熱力學、聲學等等）l矩量法 (Method of Moments) 適合于細線

2、、平面形狀結構的電磁場問題Beihang University電磁學基本方程l麥克斯韋方程組( Maxwell Equations)James Clerk Maxwell (18311879)(法拉第感應定律)(安培環(huán)路定律)(高斯定律電場)(高斯定律磁場)l物質(zhì)本構關系 (Constitutive Relations)（歐姆定律）（電極化）（磁化）JE 一切電磁場和電磁波問題均可由以上方程，以及各類具體的邊界條件所決定！Beihang University麥克斯韋方程組中的運算符散度（Divergence)旋度（Curl)連續(xù)函數(shù)的偏微分運算Beihang University設有一連續(xù)函數(shù)

3、 , 現(xiàn)欲求。FDTD的基本思想( )f x( )f x( )fx連續(xù)偏微分的有限階近似時域和空間域的離散化()()( )2 f xxf xxfxx二階中心差分近似表達式：當越小時，上式的近似程度越高。x2()()( )( ).26 f xxf xxxfxfxx實際上：Beihang UniversityFDTD空間域的離散化(1)空間域的分割離散化E Ex x分量的空分量的空間間離散分布離散分布圖圖YeeYee元胞（元胞（ x, x, y, y, z)z)H Hx x分量的空分量的空間間離散分布離散分布圖圖節(jié)節(jié)點點Beihang UniversityFDTD空間域的離散化1()yzxzEHE

4、tyx YEE YEE 元胞元胞各各電電磁磁場場分量在元胞中的位置分量在元胞中的位置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.例如：例如：H Hz z為為ExEx和和EyEy所所環(huán)繞環(huán)繞。Beihang UniversityFDTD空間偏微分的近似1()yzxzEHEtyx ( ,1)( , )

5、2()2xxxEEi jEi jyy(1, )( , )2()2yyyEEijEi jxx以以HzHz為為例：例：類似地，可實現(xiàn)各電磁場分量的空間偏微分計算。Beihang UniversityFDTD時間偏微分的近似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時間時間上的推移上的推移(1)(1)電場電場在在時間時間上取整上取整數(shù)數(shù)倍倍的的 t; t; t=n *t;(2)(2)磁磁場場在在時間時間上取上取( (整整數(shù)數(shù)1/2)1/2)倍的倍的 t; t; t=(n +1/2)*

6、t;1/21/21()2()2nnnnyzzxzEHHEtyx Beihang University麥克斯韋方程的離散化近似1()yzxzEHEtyx 以以HzHz為為例：例：1/21/2( , )( , )(1, )( , )1( ,1)( , )() nnnnnnzzyyxxzHi jHi jtEijEi jEi jEi jyx 1/21/2( , )( , )(1, )( , )( ,1)( , )() nnnnnnzzyyxxzHi jHi jEijEi jt Ei jEi jyx 上式即為Hz的更新方程，由前一時刻的磁場和前半時刻的臨近空間格點的電場即可求出最新時刻的磁場。Beiha

7、ng University麥克斯韋方程的離散化近似采取采取類類似的步似的步驟驟，可以推，可以推導導出其出其它場它場量的更新表量的更新表達達式：式：11/21/21/21/2( , )( , )( , )(1, )( , )( ,1)() zzynnnnyxxznnEi jEi jHi jHijtHi jHi jxy 例如，例如，對對于于EzEz：1()yzxzHEHtxy Beihang UniversityFDTD的離散參數(shù)的選擇元胞尺寸：元胞尺寸：邊長邊長小于最短波小于最短波長長的的1/10, 1/10, 以以減減小小數(shù)值數(shù)值色散。色散。maxmaxmaxmin1, and 10 xyz

8、222022, ,sin ()sin ()22()()22x y zktct數(shù)值色散方程：222220()xyzckkk理想色散方程：要求：要求：02k02tmax1 10 tf例如，取例如，取例如，取例如，取Beihang UniversityFDTD的離散參數(shù)的穩(wěn)定性條件 時間時間步步長長：Courant Courant 穩(wěn)穩(wěn)定性定性條條件件2221111 tcxyz對對于非色散介于非色散介質(zhì)質(zhì)，時間時間步步長長不能大于以下表不能大于以下表達達式：式：21/21/222, ,04sin (/2)()( )0,rx y zkZZZct（von Neumann method ）為為了保持了保持

9、穩(wěn)穩(wěn)定性，定性，該該方程的所有解的模必方程的所有解的模必須須小于小于1 1。域域數(shù)值數(shù)值色散方程：色散方程：Beihang University介質(zhì)電磁參量的設定( , )( , )1()yzxzHEHxi ji jty ( , )( , )1()yzxzi jiyjEHEtx 不同的元胞的電磁參量應設置為所在空間所代表的介質(zhì)的介電常數(shù)和磁導率。長方體，物質(zhì)I長方體，物質(zhì)II球體，物質(zhì)III空氣場量與介質(zhì)參數(shù)要對應場量與介質(zhì)參數(shù)要對應Beihang University色散介質(zhì)的FDTD模擬202200()( )2srj 以以Lorentz 介質(zhì)為例：介質(zhì)為例：1220120121220120

10、12ee() ,eejjrjjaa Za Zaaabb Zb ZbbbDE1201210012()nnnnnnbbbaaaDDDEEE更新方程：更新方程：近似，近似，求系數(shù)求系數(shù)Z-變換0()()()rDE0( )( )* ( )rtttDE02220002200()( )2( )sin() ( )tsrttet U t 時時域：域：頻頻域：域：Beihang University色散介質(zhì)的FDTD模擬模擬模擬Lorentz色散介質(zhì)的不同方法：色散介質(zhì)的不同方法：012012,a a a b b bMSE approach :( )( )rrError 12010012120()/ ()/nn

11、nnnnbbbaaaEDDDEEBeihang UniversityFDTD編程流程主循主循環(huán)環(huán)初始化初始化輸輸出出結結果果Beihang University編程舉例1：一維FDTD問題x x基本旋度方程：基本旋度方程：X向電導率X向磁電導率X向電流X向磁流Beihang University編程舉例1：一維FDTD問題Matlab程序代碼：% Define initial constantseps_0 = 8.854187817e-12; % permittivity of free space mu_0 = 4*pi*1e-7; % permeability of free space

12、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(domain_size/dx); % numbe

13、r of cells in 1D problem spacesource_position = 0.5; % position of the current source Jz 1.定義物理常量2.定義問題的參量和結構尺寸Beihang University% 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+1,1);eps_r_z = ones (nx+1,1); % free

14、 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 編程舉例1：一維FDTD問題（續(xù)）3.初始化場量和介質(zhì)參量陣列電場電電場電流部分流部分磁磁場場磁流部分磁流部分Beihang University% Calculate FDTD

15、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_e_z); Chyh = (2 * mu_r_y * mu_0 - dt * sigma_m_y) . ./(2 * mu_r_y *

16、 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); 編程舉例1：一維FDTD問題（續(xù)）4.計算更新方程系數(shù)電場電場部分部分磁場部分Beihang University% Define the Gaussian source waveform time = dt*0:number_of_time_steps-1.;Jz_waveform = exp(-(t

17、ime-2e-10)/5e-11).2);source_position_index = round(nx*source_position/domain_size)+1; % Subroutine to initialize plotting initialize_plotting_parameters; 編程舉例1：一維FDTD問題（續(xù)）5.定義場源源波形為高斯型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

18、 -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_positions,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,

19、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,center,fontweight,bold);text(1,1,1.1,PEC,horizontalalignment,center,fontweight,bold);Beihang University% FDTD loopfor time_step = 1:numbe

20、r_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) . + 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) . + Cez

21、hy(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 = 1 m % Subroutine to plot the current state of the fields plot_fields;end編程舉例1：一維FDTD問題（續(xù)）7.主循環(huán)程序引入電流源更新磁場更新電場設置完美電導層跟新圖面Beihang Univ

22、ersity% 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_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);dr

23、awnow;編程舉例1：一維FDTD問題（續(xù)）7.更新畫面運運行程序行程序: :fdtd_1d_code.mfdtd_1d_code.mBeihang University開放模擬空間的模擬完美匹配層（完美匹配層（Perfectly Matched Layer)的概念的概念吸波材料（損耗介質(zhì)）Beihang University完美匹配層（完美匹配層（Perfectly Matched Layer)材料的特點材料的特點*FmFmDmj0*FmFmHmj0, ,mx y zBeihang University形成 PML必須符合的兩個條件:0m*Fx*Fx 從計算空間到PML層的阻抗必須一直匹配

24、:2. 某方向的電磁參量必須是另一個方向電磁參量的倒數(shù)*Fx1*Fy*Fx1*Fy完美匹配完美匹配層層（Perfectly Matched Layer)Perfectly Matched Layer)材料的特點材料的特點Beihang University下面的材料就符合這些條件：FmFm1Dm0Hm0D00m*Fx*Fx1(x) / j01(x) / j01阻抗始終是相等的：電導率隨著PML的層數(shù)號（x)而逐漸增加,提高吸收外行電磁波的能力。完美匹配完美匹配層層（Perfectly Matched Layer)Perfectly Matched Layer)材料的特點材料的特點Beihang

25、 University完美匹配完美匹配層層PMLPML的的驗證實驗證實例例中心有一光源，四周被中心有一光源，四周被PMLPML環(huán)繞環(huán)繞可以看到外形可以看到外形電電磁波被磁波被完美吸收，完美吸收，沒沒有反射。有反射。運運行程序行程序: :fd2d_pml_code.mfd2d_pml_code.mBeihang UniversityFDTDFDTD的的應應用用實實例例(1)電磁波在損耗介質(zhì)中的傳播(2)電磁波在色散介質(zhì)中的傳播損耗介質(zhì)：4r0.0004 /S m運運行程序行程序: :FD1D_Loss.mFD1D_Loss.m只有幅值的減小，沒有色散。只有幅值的減小，沒有色散。同時有幅值的減小，以及色散。同時有幅值的減小，以及色散。020.0001 /S mDebye色散介質(zhì)0.591 10 s 01rjj運運行程序行程序: :FD1D_Dispersion.mFD1D_Dispersion.mBeihang