下载文档原格式
时域有限差分法对平面TE波的MATLAB仿真姓名：王云璐学号：2011302021（西北工业大学电子信息学院08041103，陕西西安，710072）摘要：本文分析FDTD算法的基本原理及两种典型边界条件的算法特点，通过MATLAB程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。
得到了相应的磁场幅值效果图。
最后得出用Matlab语言对FDTD算法编程的几点结论。
关键词：FDTD；MATLAB；PML；平面TE波1引言电磁场的有限差分一般是在时域上进行的，随着计算机技术的发展和应用，近年来，时域计算方法越来越受到重视。
时域有限差分法具有简单、结果直观、网格剖分简单等优点。
近些年FDTD发展的十分迅速，在许多方面都有重要应用，包括天线设计，微波电路设计，电磁兼容分析，电磁散射计算，光子学应用等等。
时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理，用来解决电磁脉冲在电磁介质中传播和反射问题的算法。
基本思想是：FDTD计算域空间节点采用Yee元胞的方法，同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区，这两个区域是以连接边界相连接，最外边是采用特殊的吸收边界，同时在这两个边界之间有个输出边界，用于近、远场转换；在连接边界上采用连接边界条件加入入射波，从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件，尽量消除反射波在吸收边界上的非物理性反射波。
本文主要简述了FDTD算法的基本原理，解的稳定性以及边界条件特点，并用Matlab语言进行对平面TE波进行了编程计算。
2FDTD基本原理时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化，用差分方程代替一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。
FDTD空间网格单元上电场和磁场各分量的分布如图1所示。
图1.Yee 氏网格及其电磁场分量分布电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元每个面的中心，每个磁场(电场)分量都有4个电场(磁场)分量环绕。
时域有限差分法
时域有限差分法的基本思想是用中心差商代替场量对时间和空间的一阶偏微商, 通过在时域的递推模拟波的传播过程, 从而得出场分布。
它最早由K.S.Yee 于1966 年提出，在此之后的20 年内，其研究进展缓慢，只是在电磁散射、电磁兼容领域有一些初步的应用。
自80 年代末，时域有限差分法成为电磁场数值计算的重要方法之一。
在声学数值计算中，时域有限差分法已应用于水声学、噪声控制及室内声学等方面的数值模拟。
时域有限差分法(Finite difference time domainmethod，FDTD)直接离散时域波动方程，不需要任何形式的导出方程，故不会因为数学模型而限制其应用范围。
它的差分格式中包含有介质的参量，只须赋予各网格相应的参量，就能模拟各种复杂的结构，这是时域有限差分法的一个突出优点。
另外，由于时域有限差分法采用步进法进行计算，故能很容易地实现各种复杂时域宽带信号的模拟，而且可以非常方便地获得空间某一点的时域信号波形。
时域有限差分法对平面TE波的MATLAB仿真摘要时域有限差分法是由有限差分法发展出来的数值计算方法。
自1966年Yee 在其论文中首次提出时域有限差分以来，时域有限差分法在电磁研究领域得到了广泛的应用。
主要有分析辐射条线、微波器件和导行波结构的研究、散射和雷达截面计算、分析周期结构、电子封装和电磁兼容的分析、核电磁脉冲的传播和散射以及在地面的反射及对电缆传输线的干扰、微光学元器件中光的传播和衍射特性等等。
由于电磁场是以场的形态存在的物质，具有独特的研究方法，采取重叠的研究方法是其重要的特点，即只有理论分析、测量、计算机模拟的结果相互佐证，才可以认为是获得了正确可信的结论。
时域有限差分法就是实现直接对电磁工程问题进行计算机模拟的基本方法。
在近年的研究电磁问题中，许多学者对时域脉冲源的传播和响应进行了大量的研究，主要是描述物体在瞬态电磁源作用下的理论。
另外，对于物体的电特性，理论上具有几乎所有的频率成分，但实际上，只有有限的频带内的频率成分在区主要作用。
文中主要谈到了关于高斯制下完全匹配层的差分公式的问题，通过MATLAB 程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。
得到了相应的磁场幅值效果图。
关键词：时域有限差分完全匹配层MATLAB 磁场幅值效果图目录摘要 (1)目录 (3)第一章绪论 (4)1.1 课题背景与意义 (4)1.2 时域有限差分法的发展与应用 (4)2.1 Maxwell方程和Yee氏算法 (7)2.2 FDTD的基本差分方程 (9)2.3 时域有限差分法相关技术 (11)2.3.1 数值稳定性问题 (11)2.3.2 数值色散 (12)2.3.3 离散网格的确定 (13)2.4 吸收边界条件 (13)2.4.1 一阶和二阶近似吸收边界条件 (14)2.4.2 二维棱边及角顶点的处理 (17)2.4.3 完全匹配层 (19)2.5 FDTD计算所需时间步的估计 (23)第三章MATLAB的仿真的程序及模拟 (25)3.1 MATLAB程序及相应说明 (25)3.2 出图及结果 (28)3.2.1程序部分 (28)3.2.2 所出的效果图 (29)第四章结论 (31)参考文献 (32)第一章绪论1.1 课题背景与意义20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。
时域有限差分法（FDTD 算法）时域有限差分法是1966年K.S.Yee 发表在AP 上的一篇论文建立起来的，后被称为Yee 网格空间离散方式。
这种方法通过将Maxwell 旋度方程转化为有限差分式而直接在时域求解, 通过建立时间离散的递进序列, 在相互交织的网格空间中交替计算电场和磁场。
FDTD 算法的基本思想是把带时间变量的Maxwell 旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应。
需要考虑的三点是差分格式、解的稳定性、吸收边界条件。
有限差分通常采用的步骤是：采用一定的网格划分方式离散化场域；对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；结合选定的代数方程组的解法，编制程序，求边值问题的数值解。
1.FDTD 的基本原理FDTD 方法由Maxwell 旋度方程的微分形式出发，利用二阶精度的中心差分近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上对连续电磁场数据的抽样压缩。
Maxwell 方程的旋度方程组为：E E HtHH Emt（1）在直角坐标系中，（1）式可化为如下六个标量方程：zz x y y y z xx x yzE tE yH xH E t E x H z HE t E zHy H ，zmz x y y m y z x x mx y z H tH yE xE H t H x E z E H t H z E y E （2）上面的六个偏微分方程是FDTD 算法的基础。
Yee 首先在空间上建立矩形差分网格，在时刻t n 时刻，F(x,y,z)可以写成),,(),,,(),,,(k j i F t n z k y j x i F t z y x F n（3）用中心差分取二阶精度：对空间离散：2),,21(),,21(),,,(xO xk j iF k j iF x t z y x F nnxi x 2),21,(),21,(),,,(yO yk j i F k ji F yt z y x F nnyj y 2)21,,()21,,(),,,(zO zk j i F k j i F zt z y x F nnzk z对时间离散：22121),,(),,(),,,(tO tk j i Fk j i Ftt z y x F n n tn t （4）Yee 把空间任一网格上的E 和H 的六个分量，如下图放置:oyxzEyHzExEzHxEyEyEzEx HyEzEx 图1 Yee 氏网格及其电磁场分量分布在FDTD 中，空间上连续分布的电磁场物理量离散的空间排布如图所示。
时域有限差分法的Matlab仿真
张通;孙晶
【期刊名称】《科技视界》
【年(卷),期】2017(000)003
【摘要】文章介绍了时域有限差分法的基本原理,利用matlab仿真,实现了用时域有限差分程序来计算二维问题空间中的电场分布.
【总页数】3页(P18-19,4)
【作者】张通;孙晶
【作者单位】吉首大学物理与机电工程学院,湖南吉首416000;吉首大学物理与机电工程学院,湖南吉首416000
【正文语种】中文
【相关文献】
1.基于高阶时域有限差分法与改进节点分析法混合求解复杂传输线网络瞬态响应[J], 王为;覃宇建;刘培国;周东明
2.时域有限差分法的Matlab仿真 [J], 郭春波
3.时域有限差分法和基于频域的有限元法模拟仿真在通r信车的天线布局设计的应用 [J], 李宏
4.时域有限差分法的Matlab仿真 [J], 黄明红;汪清泉
5.一种简化的分析光纤布拉格光栅的时域有限差分束传输法 [J], 戴劲草;王云明;张明德;孙小菡
因版权原因，仅展示原文概要，查看原文内容请购买。
一维有限时域差分方法(FDTD)计算中参数的选择
陈义万;李文兵;杜海霞;彭波勇
【期刊名称】《教育教学论坛》
【年(卷),期】2013(000)003
【摘要】在用一维时域有限差分方法计算波的传播时,为了防止出现发散的结果,要求时间步长△t＜Δx/v,但具体取多大的值,没有定论.通过编程试验,时间步长取△t＜0.75Δx/c,而时间步数Nt大于波从波源传到观察点所需要的时间,又小于反射波到达观察点的时间.而且计算的时间格点和位置格点全部都可以取整数.
【总页数】2页(P155-156)
【作者】陈义万;李文兵;杜海霞;彭波勇
【作者单位】湖北工业大学理学院,湖北武汉430068;湖北工业大学理学院,湖北武汉430068;第二炮兵工程大学,陕西西安710025;第二炮兵工程大学,陕西西安710025
【正文语种】中文
【中图分类】G642.4
【相关文献】
1.时域有限差分方法分析微带天线散射参数和方向图
2.快速时域有限差分方法计算矩形缺陷接地结构
3.区域分解时域有限差分方法(DD-FDTD)及其在散射问题中的应用
4.一种节约内存的局部一维减缩时域有限差分方法
5.利用时域有限差分方法计算雷电水平电场
因版权原因，仅展示原文概要，查看原文内容请购买。
时域有限差分法仿真一维TE 波在分裂场完全匹配层吸收边界条件下的传输一、时域有限差分法 (FDTD, Finite-Difference Time-Domain)FDTD 是1966年K.S.Yee 发表在AP 上的一篇论文建立起来的，后被称为Yee 网格空间离散方式核心思想是把带时间变量的Maxwell 旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应号称目前计算电磁学界最受关注，最时髦的算法，但还在发展完善之中国外已有多种基于FDTD 算法的电磁场计算的软件:XFDTD 等。
二、差分算法x h ∆= ()()f f x h f x ∆=+-0()()()()lim ()()()()=2x df f x f x f x h f x dx x x hf x f x h hf x h f x h h∆→∆∆+-=≈=∆∆--=+--前向差分 后向差分中心差分 222()1()()()2!df x d f x f x h f x h h dx dx +=+++222()1()()()2!df x d f x f x h f x h h dx dx -=-++333()2()()()23!df x d f x f x h f x h h h dx dx +--=++三、时域有限差分算法步骤(1)采用一定的网格划分方式离散化场域；(2)对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；(3)结合选定的代数方程组的解法，编制程序，求边值问题的数值解。
四、吸收边界条件1、问题的提出在电磁场的辐射和散射问题中，边界总是开放的，电磁场占据无限大空间，而计算机内存是有限的，所以只能模拟有限空间。
即：时域有限差分网格将在某处被截断。
这要求在网格截断处不能引起波的明显反射，因而对向外传播的波而言，就像在无限大的空间传播一样，一种行之有效的方法是在截断处设置一种吸收边界条件。
使传播到截断出的波被边界吸收而不产生反射。
2、良匹配层基本原理1994年Beernger 提出一种新颖的由吸收媒质构成的良匹配层(PML)的概念,这种人工设计的良匹配层由有耗,导电,导磁媒质组成,可吸收任意入射角,任意频率,任意偏振态的入射电磁波,且透射波以原来波速传播,且有相同的特征阻抗,在垂直入射方向衰减。
根据推导，只要这种媒质的结构参数满足一定条件,这种媒质就可以起到完全吸收入射波的作用。
五、计算所用到的电磁场公式根据麦克斯韦方程组和时域差分法，得到以下公式：()()()()1/21/21/21/2nn n n y y x x H k H k E k E k t x +-+---=∆∆()()()()11/21/201/21/211n n n n y y x x H k H k E k E k txμ++++--+-=-⋅∆∆12=1/21/21/2()()[(1/2)(1/2)]n n n nx x y y rE k E k H k H k ε++=++--1/21/21(1/2)=(1/2)-[(1)()]2n nn n y y x x H k H k E k E k +++++-六、仿真结果timeE -f i e l di n p u tfrequencyp o i n t s /w a v e l e n g t h-4time [ns]R e f l e c t e d f i e l dfrequency [GHz]R e f l e c t i o n i n d Btime [ns]T r a n s m i t t e d f i e l d-3frequency [GHz]T r a n s m i s s i o n i n d B七、源代码%*********************************************************************** % 时域有限差分法仿真一维TE波在分裂场完全匹配层吸收边界条件下的传输%*********************************************************************** %%% Eref Ein Erans% PEC PML 介质PML PEC% +-----+---+------+---------------------+-----+------+-->% z=0 z=Lz% r=1 Nref Nin Ntrans r=Nz+1%%*********************************************************************** % 物理常量mu0 = 4e-7 * pi; % 磁导率c0 = 299792458; % 光速eps0 = 1/c0^2/mu0; % 电导率eta0 = sqrt(mu0/eps0); % 波导率GHz = 10^9;mm = 10^(-3);%*********************************************************************** % 参数初始化%*********************************************************************** Lz = 1; % 长度单位：米Tmax = 4*Lz/c0; % 时间最大值Dz = 1*mm; % 空间网格尺寸R = 1/sqrt(3); %每一单位时间计算的网格单元数freq = 5*GHz; % 激励源的中心频率fmax = 9*GHz; % 画图的最大频率fmin = 1*GHz;Nz = round(Lz/Dz); % z轴上的单元格数量（round 四舍五入）n_pict = round(Nz/20); % 每n_pict步画一次图linew = 2; % 画图的线宽Dt = Dz*R/c0; % 时间步长Nt = round(Tmax/Dt); % 时间步长数量Nfft = max(2^(round(log2(Nt)+1)),2*1024); %做傅里叶变换点的数量z = linspace(0,Lz,Nz+1); % Z轴坐标t = Dt*[0:Nt-1]'; % 时间lambda = c0/freq; % 中心波长的激励源ppw = lambda/Dz % 中心频率处每一段波长的采样点omega = 2.0*pi*freq; % 角频率Nabc = 1*10; % PML边界左边和右边的层数Nin = 4; % 表面散射场Nref = 2; % 由近到远场表面Ntrans = Nz-2;Labc = Nabc*Dz; % PML边界长度%*********************************************************************** % 材料参数%*********************************************************************** Nmedia=3; % 材料数量media_par = [1.0 0.0; % 真空：电容率和电导率3.4 1*10^(-12); %3.0 0];objects = 0; % 介质的数量obj_z = Dz/2+[0.500 0.504; % Z1和Z2之间的介质10.624 0.628;0.24 0.26];obj_m = [2 2 1]; % 介质的材料obj_c = 'rrb'; % 介质的颜色epsr = ones(Nz-1,1);sig = zeros(Nz-1,1);obj_ind = [];obj_ind(:,1) = ceil(obj_z(:,1)/Dz);obj_ind(:,2) = floor(obj_z(:,2)/Dz);obj_frac(:,1) = obj_ind(:,1)-obj_z(:,1)/Dz;obj_frac(:,2) = -obj_ind(:,2)+obj_z(:,2)/Dz;for n=1:objectsepsr(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),1);sig(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),2);end%*********************************************************************** % 分配场矢量%*********************************************************************** Ex = zeros(Nz+1 ,1);Hy = zeros(Nz,1);Eref = zeros(Nt,1);Etrans = zeros(Nt,1);% 吸收边界条件E1abc = zeros(Nabc+1,1); % 左端E2abc = zeros(Nabc+1,1); % 右端H1abc = zeros(Nabc,1);H2abc = zeros(Nabc,1);zpmlH = linspace(0,Labc,Nabc)'/Labc;zpmlE = linspace(Dz/2,Labc-Dz/2,Nabc-1)'/Labc;abc_pow = 3;sigma_max = -(abc_pow+1)*log(10^(-4))/(2*eta0*Labc);sigmaH2 = sigma_max*zpmlH.^abc_pow*mu0/eps0;sigmaE2 = sigma_max*zpmlE.^abc_pow;sigmaE1 = sigmaE2((Nabc-1):-1:1);sigmaH1 = sigmaH2(Nabc:-1:1);figure(4)plot(sigmaH2)%***********************************************************************% 波激励%***********************************************************************Hdelay = Dt/2-Dz/2/c0;% 从左侧加入射脉冲tau = 80.0e-12;delay = 2.1*tau;Ein = sin(omega*(t-delay)).*exp(-((t-delay).^2/tau^2));Hin = 1/eta0*sin(omega*(t-delay+Hdelay)).*exp(-((t-delay+Hdelay).^2/tau^2));% 谐波从左侧入射的时间Einmax = max(abs(Ein));%***********************************************************************% 时间步进%***********************************************************************for n = 1:Nt;% 从-Dt/2到Dt/2更新磁场Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1); % 主网格Hy(Nin) = Hy(Nin) - (Dt/mu0/Dz) * Ein(n);H1abc = (H1abc.*(mu0/Dt-sigmaH1/2) - diff(E1abc)/Dz) ./ (mu0/Dt+sigmaH1/2); % 从ABC向左H2abc = (H2abc.*(mu0/Dt-sigmaH2/2) - diff(E2abc)/Dz) ./ (mu0/Dt+sigmaH2/2); % 从ABC向右% 从0到Dt更新EEx(2:Nz) = (Ex(2:Nz).*(eps0*epsr/Dt-sig/2) - diff(Hy,1)/Dz)./(eps0*epsr/Dt+sig/2); % 除去边界的主网格Ex(Nin) = Ex(Nin) - Hin(n)/(Dz*eps0/Dt);Ex(1) = Ex(1) - (Dt /eps0) * (Hy(1)-H1abc(Nabc))/Dz; %左边ABC和主网格之间的边界E1abc(Nabc+1) = Ex(1); % 补充一点，以简化方案E1abc(2:Nabc) = (E1abc(2:Nabc).*(eps0/Dt-sigmaE1/2) - diff(H1abc,1)/Dz) ./ (eps0/Dt+sigmaE1/2);Ex(Nz+1) = Ex(Nz+1) - (Dt /eps0) * (H2abc(1)-Hy(Nz))/Dz;E2abc(1) = Ex(Nz+1);E2abc(2:Nabc) = (E2abc(2:Nabc).*(eps0/Dt-sigmaE2/2) - diff(H2abc,1)/Dz) ./ (eps0/Dt+sigmaE2/2);% 在选定点采样电场if mod(n,n_pict) == 0figure(1);subplot(2,1,1);plot(z,Ex,'LineWidth',1);title(['time is ',num2str(n*Dt*10^9),'ns'])axis tight;ax=axis;axis([ax(1:2) -1.2*Einmax 1.2*Einmax]);ax=axis;ylabel('E-field');for p=1:objectsline([obj_z(p,1) obj_z(p,1)], ax(3:4),'color',obj_c(p));line([obj_z(p,2) obj_z(p,2)], ax(3:4),'color',obj_c(p));line([obj_z(p,1) obj_z(p,2)], [ax(3) ax(3)],'color',obj_c(p),'LineWidth',4);line([obj_z(p,1) obj_z(p,2)], [ax(4) ax(4)],'color',obj_c(p),'LineWidth',4);endsubplot(2,1,2);plot([E1abc(:)' Ex' E2abc(:)'],'LineWidth',1)title(['time is ',num2str(n*Dt*10^9),'ns'])ax = axis;line([Nabc Nabc], ax(3:4),'color','r');line(Nin+[Nabc Nabc], ax(3:4),'color','b');line(Ntrans+2+[Nabc Nabc], ax(3:4),'color','b');line([Nz+2+Nabc Nz+2+Nabc], ax(3:4),'color','r');axis tight;ylabel('E-field');%plot(Eabc(2:Nabc,1:2));pause(0.01);endEref(n) = Ex(Nref);Etrans(n) = Ex(Ntrans);end%***********************************************************************% 后处理%***********************************************************************Ehin = fft(Ein,Nfft)/Nt*2;Df = 1/(Dt*Nfft);freqN = round(fmax/Df);freqN1 = round(fmin/Df);f = linspace(Df*freqN1,Df*freqN,freqN-freqN1+1)/10^9;t = t*10^9; % 纳米秒（ns）figure(2);% 入射波clf;subplot(2,3,1);plot(t,Ein,'r','LineWidth',linew)axis tight;xlabel('time');ylabel('E-field');subplot(2,3,4);warning off;A=plotyy(f,abs(Ehin(freqN1:freqN)),f,c0./f/Dz/10^9);warning on;set(A(2),'yLim',[0 40],'yTick',[0:5:40],'xlim',[fmin fmax]/10^9,'xGrid','on','yGrid','on','LineWidth',linew); set(A(1),'xlim',[fmin fmax]/10^9,'LineWidth',linew);set(get(A(2),'Ylabel'),'String','points/wavelength');set(get(A(1),'Ylabel'),'String','input');set(get(A(1),'Xlabel'),'String','frequency');% 反射subplot(2,3,2);plot(t,Eref,'LineWidth',linew); axis tight;xlabel('time [ns]');ylabel('Reflected field');subplot(2,3,5);Ehref = fft(Eref,Nfft)/Nt*2;plot(f,20*log10(abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;ax=axis;axis([fmin/10^9 fmax/10^9 ax(3:4)]);axis tight;xlabel('frequency [GHz]');ylabel('Reflection in dB');% 传播subplot(2,3,3);plot(t,Etrans,'LineWidth',linew); axis tight;xlabel('time [ns]');ylabel('Transmitted field');subplot(2,3,6);Ehtrans = fft(Etrans,Nfft)/Nt*2;plot(f,20*log10(abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;ax=axis;axis([fmin/10^9 fmax/10^9 ax(3:4)]);xlabel('frequency [GHz]');ylabel('Transmission in dB');tr=interp1(f,abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN)),5) re=interp1(f,abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN)),5) tr^2+re^2。
