时域有限差分方法林志立实用教案

1、电磁学基本(jbn)方程 麦克斯韦(mi k s wi)方程组( Maxwell Equations)James Clerk Maxwell (18311879)(法拉第感应(gnyng)定律)(安培环路定律)(高斯定律电场)(高斯定律磁场)l物质本构关系 (Constitutive Relations)（欧姆定律）（电极化）（磁化）JE 一切电磁场和电磁波问题均可由以上方程，以及各类具体的边界条件所决定！第1页/共28页第一页，共28页。麦克斯韦(mi k s wi)方程组中的运算符散度（Divergence)旋度（Curl)连续函数的偏微分运算(yn sun)第2页/共28页第二页，共28

2、页。设有一连续函数 , 现欲求(y qi)。FDTD的基本(jbn)思想( )f x( )f x( )fx连续偏微分的有限(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、空间间离散分布离散分布图图节节点点第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所所环绕环绕。第5页/共28页第五页，共28页。FDTD空间(kngjin)偏微分的近似1()yzxzEHEtyx ( ,1)( , )2()2xxxEEi jEi jyy(1, )( , )2()2yyyEEijEi jxx以HzHz为例：类似地，可实现各电磁场分量(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)在在时间时间上上取整取整数数倍的倍的t; t; t=n t=n * *t;t;(2)(2)磁磁场场在在时间时间上取上取( (整整数数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的离散参数(cnsh)的选择元胞尺寸元胞尺寸(ch cun)：边长小于最短波长的：边长小于最短波长的1/10, 以减小数值色散。以减小数值色散。maxmaxmaxmin1, and 10 xyz 222022, ,sin ()sin ()22()()22x y zktct数值(shz)色散方程：222220()xyzckkk理想色散方程：要求：要求：02k02tmax1 10 tf例如，取例如，取例如，取例如，取第10页/共28页第十页，共28

8、页。FDTD的离散(lsn)参数的稳定性条件 时间时间(shjin)(shjin)步步长长：Courant Courant 稳稳定性定性条条件件2221111 tcxyz对对于非色散于非色散(ssn)(ssn)介介质质，时间时间步步长长不能大于以下表不能大于以下表达达式：式：21/21/222, ,04sin (/2)()( )0,rx y zkZZZct（von Neumann method ）为为了保持了保持稳稳定性，定性，该该方程的所有解的模必方程的所有解的模必须须小于小于1 1。域数值色散方程：第11页/共28页第十一页，共28页。介质电磁(dinc)参量的设定( , )( , )1(

9、)yzxzHEHxi ji jty ( , )( , )1()yzxzi jiyjEHEtx 不同的元胞的电磁(dinc)参量应设置为所在空间所代表的介质的介电常数和磁导率。长方体，物质I长方体，物质II球体，物质III空气场量与介质参数场量与介质参数(cnsh)要对应要对应第12页/共28页第十二页，共28页。色散(ssn)介质的FDTD模拟202200()( )2srj 以以Lorentz 介质介质(jizh)为例：为例：122012012122012012ee() ,eejjrjjaa Za Zaaabb Zb ZbbbDE1201210012()nnnnnnbbbaaaDDDEEE更新

10、更新(gngxn)方程：方程：近似，近似，求系数求系数Z-变换0()()()rDE0( )( )* ( )rtttDE02220002200()( )2( )sin() ( )tsrttet U t 时时域：域：频频域：域：第13页/共28页第十三页，共28页。色散介质(jizh)的FDTD模拟模拟模拟Lorentz色散色散(ssn)介质的不同方法：介质的不同方法：012012,a a a b b bMSE approach :( )( )rrError 12010012120()/ ()/nnnnnnbbbaaaEDDDEE第14页/共28页第十四页，共28页。FDTD编程流程(lichng

11、)主循主循环环(xnhun)(xnhun)初始化初始化输输出出(shch)(shch)结结果果第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.定义问题的参量和结构(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问题（续）问题（续）3.初始化场量和介质(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问题（续）问题（续）4.计算更新方程(fngchng)系数电场电场(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问题（续）问题（续）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问题（续）7.主循环(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问题（续）7.更新(gngxn)画面运运行行(ynxng)(ynxng)程序程序: :fdtd_1d_code.mfdtd_1d_code.m第22页/共28页第二十二页，共28页。开放(kifng)模拟空间的模拟完美匹配完美匹配(ppi)层（层（Perfectly Matched Layer)的概念的