科技1101班-杨棣.doc

开放式光电子器件仿真实验 基于FDTD一维电场、磁场设计 班 级: 科技1101姓 名： 杨棣成绩指导教师一、实验目的和选题背景（一）、实验目的1. 掌握麦克斯韦方程组及其fdtd形式。2. 掌握麦克斯韦方程组的差分离散方法从而进行差分离散得到一维的时域简单分析。3. 掌握FDTD中Yee元胞的电场和磁场分析。（二）、选题背景FDTD时域有限差分法 是1966年K.S.Yee发表在AP上的一篇论文建立起来的，后被称为Yee网格空间离散方式, 它的核心思想是把带时间变量的Maxwell旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应。号称目前计算电磁学界最受关注，最时髦的算法，但还在发展完善之中，国外已有多种基于FDTD算法的电磁场计算的软件，但在国内基于MATLAB的FDTD算法才刚起步没多少年，而且研究的人也不多。关键的三大要素就是：差分格式、解的稳定性、吸收边界条件；FDTD的特点：广泛的应用性、节约运算和存储空间、适合并行计算、计算程序的通用性、简单直观，容易掌握。二、实验原理2.1FDTD基本原理1、设函数f(x), 独立变量x有很小的增量 ，则 （1） （2）这三种不同的近似表达式引入的误差，可由泰勒级数来分析。有Taylor公式可得: （3）Maxwell方程组即 （4）maxwell旋度方程可以推出此六个耦合方程 （5）式中，时电流密度，反映电损耗，是磁流密度，单位，反映磁损耗。主要与上式对应。各向同性介质中的本构关系： 令代表在直角坐标中的任何一个分量，离散符号取为 （6）关于时间和空间的一阶偏导数取中心差分近似为 （7）可以看出，每一节点上沿某一方向场分量的一阶偏微分可以用在该方向上相邻两点的一阶中心差商来描述，将式(1)用一阶中心差商方程取代，整理后便得到一阶差分方程，它具有二阶精度3。 2、 有限差分法通常采用的步骤：（1）采用一定的网格划分方式离散化场域（2）对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组（3）结合选定的代数方程组的解法，编制程序，求边值问题的数值解（4）FDTD就是根据这个步骤结合自身的特点进行的。2.2 FDTD的Yee元胞Yee首先在空间上建立矩形差分网格，在时刻nt时刻，F(x,y,z)可以写成： ，用中心差分取二阶精度：对空间离散:（8）对时间离散: （9）为了满足(1)式空间精度的要求，并满足(2)式，Yee把空间任一网格上的E和H的六个分量E,H场分量，取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理 图1 Yee模型如图1所示，Yee单元有以下特点2：1）E与H分量在空间交叉放置，相互垂直；每一坐标平面上的E分量四周由H分量环绕，H分量的四周由E分量环绕；场分量均与坐标轴方向一致。2）每一个Yee元胞有8个节点，12条棱边，6个面。棱边上电场分量近似相等，用棱边的中心节点表示，平面上的磁场分量近似相等，用面的中心节点表示。3）每一场分量自身相距一个空间步长，E和H相距半个空间步长4）每一场分量自身相距一个时间步长，E和H相距半个时间步长，电场取n时刻的值，磁场取n+0.5时刻的值；即：电场n时刻的值由n-1时刻的值得到，磁场n+0.5时刻的值由n-0.5时刻的值得到；电场n时刻的旋度对应n+0.5时刻的磁场值，磁场n+0.5时刻的旋度对应(n+0.5)+0.5时刻的电场值，逐步外推。5）3个空间方向上的时间步长相等， 以保证均匀介质中场量的空间变量与时间变量完全对称。应用这种离散方式，将含时间变量的Maxwell方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。由电磁问题的初值和边界条件，就可以逐步推进地求解以后各时刻空间电磁场分布。2.3一维问题 均匀平面波（TEM波）是一维问题，设电磁波沿z轴方向传播，则，场量和介质参数均与x，y无关，即，麦克斯韦方程为 （7）和 （8）旋转坐标轴后可以只保留一组公式4，设保留（7）Yee元胞如图2所示 图2 一维Yee元胞差分格式为 （9） （10） 如果介质无损耗，则一维空间的FDTD模拟的电场、磁场的计算程序设计及模拟结果% Define initial constanteps_0 = 8.854187817e-12; % permittivity of free space mu_0 = 4*pi*1e-7; % permeability 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(domain_size/dx); % number of cells in 1D problem spacesource_position = 0.5; % position of the current source Jz% 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 spacesigma_e_z = zeros(nx+1,1); % free spaceChyh = 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% 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_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);% 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;% 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) . + 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: