时域有限差分法论文
时域有限差分法中的不均匀性问题研究
ｄｉｆｆｅｒｅｎｔ ｓｉｄｅｓ ｏｆ ｔｈｅ ｉｎｔｅｒｆａｃｅ，ａｒｅ
ｉｎｖｅｓｔｉｇａｔｅｄ
ｇｒｉｄｓ，ａｎｄ
ｔｈｅ ａｃｃｕｒａｃｙ ｏｒｄｅｒ ｏｆ ｔｈｅ ａｃｒｏｓｓ－ｉｎｔｅｒｆａｃｅ
ｅｑｕａｔｉｏｎｓ
ａｌｔｅｒｎａｔｉｎｇ－
ｄｉｆｆｅｒｅｎｔ
ｅｒｒｏｒ
ｄｉｒｅｃｔｉｏｎ ｉｍｐｌｉｃｉｔ ｆｉｎｉｔｅ—ｄｉｆｆｅｒｅｎｃｅ ＡＤＩ－ＦＤＴＤ ｆｏｒｍｕｌａｓ
ａｒｅ
ｔｉｍｅ—ｄｏｍａｉｎ（ＡＤＩ—ＦＤＴＤ）ｍｅｔｈｏｄ，ｔｗｏ
ｃａｎ
ｉｎｔｒｏｄｕｃｅｄ，ｂｏｔｈ ｏｆ ｗｈｉｃｈ
ｄｅｃｒｅａｓｅ ｔｈｅ ｒｅｆｌｅｃｔｉｏｎ
Ｋｅｙｗｏｒｄ：
ＦＤＴＤ
ｉｎｈｏｍｏｇｅｎｅｏｕｓ ｍｅｄｉｕｍ
ＡＤＩ．ＦＤＴＤ
ｅｌｅｃｔｒｉｃａｌ ｕｎｉｆｏｒｍ
ｇｒｉｄｓ
ＰＭＬ
ｍｅｄｉｕｍ
独创性（或创新性）声明
本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研究 成果。尽我所知，除了文中特别加以标注和致谢中所罗列的内容以外，论文中不 包含其他人已经发表或撰写过的研究成果；也不包含为获得西安电子科技大学或 其它教育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做 的任何贡献均已在论文中做了明确的说明并表示了谢意。 申请学位论文与资料若有不实之处，本人承担一切相关责任。
ｅｑｕａｔｉｏｎｓ ｉｓ ｐｒｅｓｅｒｖｅｄ ｓｅｃｏｎｄ・－ｏｒｄｅＬ
Ｗｉｔｈ
ａｌｗａｙｓ，ｔｈｅ ｈｉｇｈｅｒ ｔｈｅ ａｃｃｕｒａｃｙ ａｔ ｔｈｅ ｍｅｄｉｕｍ ｉｎｔｅｒｆａｃｅ，ｔｈｅ ｓｍａｌｌｅｒ Ｆｏｒ ｔｈｅ
微带间隙和宽度跃变的时域有限差分法分析
微带间隙和宽度跃变的时域有限差分法分析【摘要】本文应用时域有限差分法分析微带间隙和宽度跃变两种不连续性，准确地给出了它们的散射参数，并与商用软件IE3D的计算结果做了对比，两者的计算结果基本一致。
【关键词】时域有限差分法；微带间隙；微带宽度跃变【Abstract】Microstrip gap and microstrip step-in-width is analyzed using FDTD in this article，and S-parameter is obtained. It is almost the same to the result of calculation using IE3D.【Key words】FDTD；Microstrip gap；Microstrip step-in-width0 引言在微波电路中，通常会遇到微带传输线的不连续性。
这些不连续性的出现，都给电路设计带来新的课题。
因为不连续性对电路带来的负面影响或利用不连续性的影响提高电路的性能等，而成为电路设计中必须考虑的因素。
本文主要研究微带间隙和宽度跃变两种不连续性。
对于微带线间隙和宽度跃变，已经有很多文献研究，如文献[1][2]。
根据静态分析法和全波分析法得到的结果通过曲线拟合表示成解析形式，已经被广泛应用在微带电路设计中[3-5]。
但这些经验公式多用在频率较低的情况下，频率较高时误差就比较大甚至不适用。
因此，本文选用FDTD法分析微带间隙和宽度跃变，通过一次计算就可给出丰富的宽频带信息[6]。
1 基本原理、公式时域有限差分法的原理非常简单，就是直接将时域Maxwell方程组的两个旋度方程中关于空间变量和时间变量的偏导数用差商近似，从而转换为离散网格点上的时域有限差分方程。
在各向同性，无耗条件下的麦克斯韦方程为[7]：■=■■×H■=-■■×E （1）常规区域x、y、z方向的Maxwell方程变为：？滋■=-（■-■）？着■=■-■？滋■=-（■-■）？着■=■-■ （2）？滋■=-（■-■）？着■=■-■PML区域，Maxwell方程修正为：？滋■+？滓■H■=-■ ？着■+？滓■E■=■？滋■+？滓■H■=■ ？着■+？滓■E■=-■？滋■+？滓■H■=-■ ？着■+？滓■E■=■（3）？滋■+？滓■H■=■ ？着■+？滓■E■=-■？滋■+？滓■H■=-■ ？着■+？滓■E■=■？滋■+？滓■H■=■ ？着■+？滓■E■=-■图1 FDTD离散中的Yee元胞采用图1所示的Yee的离散网格，利用二阶精度的中心差分近似把方程（2）和（3）中的公式转换为差分形式，就得到TD-FD法迭代公式，然后编程进行分析计算。
时域有限差分法二维
时域有限差分法二维1. 引言时域有限差分法（Finite Difference Time Domain, FDTD）是一种常用的数值计算方法，用于求解电磁场在时域中的传播和辐射问题。
本文将以二维情况为例，深入探讨时域有限差分法的原理和应用。
通过本文的介绍和解读，您将更全面地理解这一方法，并能够灵活应用于相关领域。
2. 时域有限差分法简介2.1 原理概述时域有限差分法是一种迭代求解偏微分方程的方法，通过将时域和空间离散化，将连续问题转化为离散问题。
在二维情况下，假设空间网格分辨率为Δx和Δy，时间步长为Δt。
根据电磁场的麦克斯韦方程组，可以利用中心差分公式进行离散化计算，得到求解方程组的更新方程。
2.2 空间离散化对于二维情况，空间离散化可以采用正交网格或非正交网格。
常见的正交网格包括方形格点、Yee网格等，而非正交网格则具有更灵活的形态。
根据需要和应用场景，选择合适的离散化方法对问题进行求解。
2.3 时间离散化时间离散化主要有显式和隐式两种方法。
显式方法将时间推进方程展开成前一时刻的电场和磁场与当前时刻的源项之间的关系，容易计算但对时间步长有限制；隐式方法则是通过迭代或矩阵计算求解当前时刻的电场和磁场。
3. 时域有限差分法的应用领域时域有限差分法广泛应用于电磁场传播和辐射问题的数值模拟中。
以下是几个典型的应用领域：3.1 辐射问题时域有限差分法可以模拟电磁波在空间中的辐射传播过程。
可以用于分析天线的辐射特性，设计无线通信系统的天线，或者分析电磁波在无线电频段的传播情况。
3.2 波导问题对于波导结构，时域有限差分法可以求解其模式、传输特性等问题。
波导结构广泛应用于光子学器件、微波器件等领域，时域有限差分法为建立数值模型和解析波导特性提供了一种有效的数值计算手段。
3.3 散射问题时域有限差分法在散射问题的数值模拟中也有重要应用。
通过模拟散射体与电磁波的相互作用过程，可以研究和分析散射体的散射特性，例如雷达散射截面的计算、微波散射问题等。
时域有限差分法关键技术及其应用研究
时域有限差分法关键技术及其应用研究时域有限差分法关键技术及其应用研究1. 引言时域有限差分法（Finite Difference Time Domain, FDTD）是一种常见的数值电磁计算方法，被广泛应用于电磁场的数值模拟和分析。
本文将介绍FDTD方法的基本原理及其一些关键技术，重点探讨其在电磁场模拟、天线研究和光学器件设计等领域的应用。
2. FDTD方法基本原理FDTD方法采用时空网格来离散求解麦克斯韦方程组，通过迭代的方式计算电磁场的时变分布。
其基本原理是利用麦克斯韦方程组的时域形式，将电场和磁场的空间导数用有限差分的形式进行近似，通过时间步进来模拟电磁场的时域行为。
FDTD方法的关键是对时空网格的离散化处理。
在时域，时间和空间被离散为等间距的格点，电磁场在格点之间通过有限差分方程进行计算，从而得到电场和磁场在每个格点的数值。
通过时间步进的迭代计算，可以模拟电磁场随时间的演化过程。
3. FDTD方法的关键技术3.1 源的建立在FDTD方法中，需要设置适当的源来激发电磁场的变化。
常见的源包括点源、平面波源和边界条件处理等。
点源是在空间某一点施加突变的电场或磁场，用于模拟电磁波的辐射和传播；平面波源是在一个平面波入射，模拟平面波在介质中的传播行为；边界条件处理则是为了模拟无限大空间中的电磁波的传播。
3.2 时间步进时间步进是FDTD方法中的一个关键技术，决定了电场和磁场的更新方式。
常用的时间步进算法有显式和隐式两种。
显式时间步进是根据已知的电场和磁场的数值，通过有限差分方程计算新的电场和磁场的值；隐式时间步进则是使用迭代或矩阵求解的方法，利用已知的旧场和新场的关系求解新场。
3.3 网格约束条件FDTD方法中需要设置一些约束条件，以满足电磁场在网格边界条件下的数值计算。
常见的约束条件有吸收边界条件和周期性边界条件。
吸收边界条件是用于吸收入射电磁波的反射波，常用的吸收边界条件有Mur吸收边界条件和PML吸收边界条件；周期性边界条件是为了模拟周期性结构或周期性辐射场景，将仿真空间分割成无限个重复的周期结构。
时域有限差分方法的研究
时域有限差分方法的研究戴国强;余震虹;高磊;李科【摘要】Finite-difference time-domain method (FDTD) has developed into a mature numerical calculation method. The paper introduces the basic properties of FDTD, the perfectly matched layer boundary condition and excitation source. By setting the perfectly matched layer boundary condition, Matlab software is adopted for the simulation of the free space scattering field to solve the electromagnetic scattering problems. The simulation results in different time step are given to realize the visualization of the whole solving process.%时域有限差分法(FDTD)已经发展成为一种成熟的数值计算方法.介绍时域有限方法的一些基本性质,完美匹配层的边界条件和激励源设置.通过对完美匹配层边界条件的设置,利用Matlab软件对自由空间散射场进行了仿真,解决电磁场散射等问题,给出了在不同时间步长下的仿真结果,实现了整个求解过程的可视化.【期刊名称】《现代电子技术》【年(卷),期】2013(036)001【总页数】4页(P140-143)【关键词】时域有限差分法;电磁场;边界条件;时间步长【作者】戴国强;余震虹;高磊;李科【作者单位】江南大学物联网工程学院,江苏无锡214122;江南大学物联网工程学院,江苏无锡214122;江南大学物联网工程学院,江苏无锡214122;江南大学物联网工程学院,江苏无锡214122【正文语种】中文【中图分类】TN911.7-340 引言麦克斯韦偏微分方程组是宏观电动力学的基本方程，它完整地描述了宏观电磁现象的基本运动规律，是一切经典电磁波理论研究和应用开发的基础。
时域有限差分算法
时域有限差分算法Finite-Difference Time-Domain (FDTD) Algorithm时域有限差分算法（Finite-Difference Time-Domain，FDTD）FDTD is a numerical technique used to solve Maxwell's equations in the time domain.时域有限差分算法是一种用于在时域中求解麦克斯韦方程的数值技术。
It discretizes the spatial and temporal domains, allowing for the simulation of electromagnetic wave propagation and interaction with complex structures.该算法将空间和时间域离散化，从而能够模拟电磁波的传播以及与复杂结构的相互作用。
The algorithm is widely used in various fields, including antenna design, microwave engineering, and electromagnetic compatibility analysis.该算法广泛应用于多个领域，包括天线设计、微波工程和电磁兼容性分析。
The main advantage of FDTD is its ability to handle arbitrary geometries and material properties, making it a powerful tool for electromagnetic modeling and simulation.时域有限差分算法的主要优势在于其能够处理任意几何形状和材料属性，使其成为电磁建模和模拟的有力工具。
However, it can be computationally demanding, especially for large-scale problems, due to the need to discretize both space andtime.然而，由于需要同时离散化空间和时间，时域有限差分算法在计算上可能要求较高，尤其是对于大规模问题。
基于交变隐式差分方向方法的时域有限差分法
基于交变隐式差分方向方法的时域有限差分法摘要本文主要针对基于交变隐式差分方向方法的时域有限差分法（Alternating Direction Implicit Finite Difference Time Domain method，简称ADI-FDTD方法）做了一定研究。
论文首先介绍了二维ADI-FDTD方法，就其数值稳定性和数值色散特性进行了研究，验算了ADI-FDTD的Mur吸收边界条件对其稳定性的影响。
关键词：基于交变隐式差分方向方法的时域有限差分法有限时域差分法数值稳定数值色散吸收边界条件一、二维ADI-FDTD 方法基本原理基于交变隐式差分方向方法的时域有限差分法也就是ADI-FDTD 方法。
传统的FDTD 方法,属于显示差分方法,因此具有显示差分方法的共同特性。
其离散格式有以下两个方面的特点：一个就是数值色散对空间离散网格的要求,空间离散网格尺寸必须为所要模拟电磁波最短波长的1/12，通常在程序中取1/20以减少数值计算带来的误差。
第二就是时间和空间离散间隔之间应当满足Courant-Friedrich-Levy(CFL)稳定性条件，或者简称Courant 稳定条件，即：222max )(1)(1)(11z y x t v ∆+∆+∆≤∆其中Vmax 为电磁波在媒质中传播的最大相速。
如果时间步长不满足上述条件，将导致FDTD 的计算发散，特别是目标较之入射波的波长有细微结构时，随着空间网格尺寸变小，为了满足稳定性条件，时间步长也相应的取小，致使总的CPU 计算时间有可能达到无法实现的地步。
为了提高FDTD 方法的计算效率和应用范围，90年代后期提出了多种与其他技术相结合的混合方法。
第六章就介绍最后一种基于交变隐式差分方向方法的时域有限差分法，也就是ADI-FDTD 方法。
与显式差分方法相反，隐式差分格式总是稳定的，其时间步长仅受数值误差的限制。
但是，隐式差分格式的缺点是需要通过矩阵求逆，或者是迭代求解大型线性方程组，计算复杂且量大。
毕业论文 毕业设计】毕业论文 外文翻译 中英文 利用时域有限差分方法演示电磁场的传播
英文参考资料原文及翻译：An Interactive Demonstration of ElectromagneticWave Propagation Using Time-Domain FiniteDifkrencesAbstract：The finite difference time-domain (FDTD) method is one of the most widely used computational methods in ing FDTD, Maxwell’s equations are solved directly in the time domain via finite differences and time stepping. The basic approach is relatively easy to understand and is an alternative to the more usual frequency-domain approaches. In order to take advantage of this, an interactive personal computer program based on FDTD has been developed. The program directly solves Maxwell’s equations via finite differences. The solution is for one dimension, corresponding to normal incidence propagation through a planar stratified medium. The program displays an electromagnetic pulse as it propagates through the medium. Most textbooks devote very little space to time-domain solutions, but instead emphasize a complex time-harmonic approach. While there are good reasons for this, a timedomain demonstration provides a much more intuitive view of wave propagation. While the program can be used as a visual “bounce diagram” movie,it is much more versatile. Since Maxwell’s equations are solved directly, the reflected and transmitted pulse amplitudes demonstrate how the reflection and transmission coefficients determine reflected and transmitted wave amplitudes. Since lossy material layers can he included, frequency dispersion can be demonstrated. If the student’s background permits, an FFT can be used to transform the time domain results to the frequency domain to obtain, for example, complex reflection and transmission coefficients for a dielectric slab for comparison to results obtained via frequency domain transmission line techniques. At the graduate level, the program can be used to demonstrate the FDTD technique as a powerful tool for solving real-world electromagnetic scattering and coupling problems.I.INTRODUCTIONCL ASSROOM demonstrations using personal com- puters (PC’s) are considered to be very useful teaching tools. Even more useful are interactive PC programs which are controlled by the student and which illustrate the phenomena being considered. Programs which illustrate pulses traveling back and forth on a transmission line, Manuscript received November 30, 1988; revised April 18, 1989. R. J. Luebbers and K. S. Kunz are with the Department of Electrical Engineering, Pennsylvania State University, University Park. PA 16802. K. A. Chamberlin is with the Department of Electrical and Computer Engineering, University of New Hampshire, Durham, NH 03824. IEEE Log Number 8930378. “bounce diagram” movies, are easy to develop and readily available. However, this type of demonstration is quite limited. They assume pulses with zero rise time which nevertheless propagate without dispersion, definitely not illustrating actualphenomena. The approach presented in this paper involves a PC demonstration program which is based on an actual solution of Maxwell’s equations. The calculations a re performed via the finite difference time-domain (FDTD) technique [ 11, which actually solves Maxwell’s differential equations via finite differences. To simplify the graphical displays, the geometry has been made one dimensional, so that the solution is actually for a plane wave normally propagating through a planar stratified medium. Each layer may have its permittivity, permeability, and conductivity specified, although to provide reasonable menu choices, the geometries which can be specified interactively are limited. Even at the undergraduate level, the FDTD formulation can be presented, and in fact for most students, it is more readily understood than the usual time-harmonic approach. This lets the students understand the theoretical basis for the graphics displayed and the menu choices available. While the program can be used as a visual “bouncediagram” movie, it is much more versatile. Since Maxwell’s equations are solved directly, the reflected and transmitted pulse amplitudes demonstrate how the reflection and transmission coefficients determine reflected and transmitted wave amplitudes. Since lossy material layers can be included, frequency dispersion can be demonstrated. If the student’s background permits, an FFT is provided to transform the time-domain results to the frequency domain to obtain, for example, complex reflection and transmission coefficients for a dielectric slab for comparison to results obtained via frequency-domain transmission line techniques. At the graduate level, the program can be used to demonstrate the FDTD techniques as a powerful tool for solving real-world electromagnetic scattering and coupling problems.11. ONE-DIMENSIONFADLT D FORMULATIONIn order to reduce the complexity of programming anddisplaying FDTD computations, we formulate the difference equations in one dimension. If we arbitrarily choose to retain Ey and H, with propagation in the x direction, then Maxwell's equations become1()y z y E H E t xσε∂∂=--∂∂ (1) 1y z E H t xμ∂∂=-∂∂ (2) Following Yee's [l] notation, we let x=i ×x ∆，t=n ×t ∆,F n (i)=F(i ⨯∆x,n ∆⨯t) Still following Yee's approach, we interleave Ey and H, in space and time, and specify E, p, and (r at discrete points (actually, layers for our 1-D formulation).After approximating the differential equations as difference equations and simplifying, we easily obtain11122()11()()()[()()]()()(()())22n n n ny yz z i t E i E i H i H i i t i x i t i εεσεσ+++∆=-+--+∆∆+∆ (3) 112211()()[(1)()]22(1/2)n n n n z z y y t H i H i E i E i x i μ+-∆+=+-+-∆+ (4) Equations (3) and (4) are readily programmed for solutionby iteration of each equation alternately.The next considerations are the stability and accuracyconstraints on Ax and A t . It must be pointed out that Ax must be much less than the minimum wavelength of interest and the minimum scatterer dimension for good accuracy. For stability, time steps must be small enough so that field values can affect only nearest neighbors during one time step internal. For our one-dimensional equations, a necessary stability criterion isc t x ∆<∆ (5) where c is the speed of light.A conservative choice is0.5c t x ∆=∆ （6）and this is used in the demonstration model. We must also consider excitation of the difference equations. One can discuss different pulse shapes relative to their frequency content. A traditional choice is the Gaussian pulse. If we let Eo be the peak amplitude, xp be the original location of the peak, and 2 w - Ax be the pulse width to 0.001 Eo amplitude, then the Gaussian pulse electric field is given by20exp[()]y p E E x x ct β=--- （7） with22ln(0.001)ln(0.001)()(2)w x wc t β==∆∆ （8） Writing (7) as a difference equation suitable for exciting a propagating wave, we have, with xp locating the initial position of the pulse peak,020/()exp[ln(0.001)()]p y i x xE i E w -∆=- （9） For propagation in the plus-x direction, the corresponding magnetic field is given as1/2200/1/4(1/2)exp[ln(0.001)()]p z i x x E H i Z w-∆++=- （10） where Zo is the impedance of free space.The final consideration before programming the above equations is absorbing pulses as they are incident at the limits of the problem space.. Significant literature exists on FDTD absorbing boundary conditions, but as our problem formulation allows only normal incidence, the problem is considerably simplified over 2-D and 3-D formulations. An adequate absorbing boundary condition is given by Mur [2] which can be reduced to one dimension. Some examples of the results which can be obtained using the one-dimensional FDTD formulation are illustrated in Figs. 1-3. The geometry consisted of a one-dimensional stratified medium with 512 layers. The medium is free space except for layers 250-309 which are some type of dielectric. Layer 310 may be specified as a perfect absorber. Each layer ( A x ) is 1.5 mm, so the dielectric slab is 9 cm thick.译文：利用时域有限差分方法演示电磁场的传播摘要——有限差分时域分析法（FDTD）是一种在电磁学中运用最广泛的计算方法。
基于时域有限差分法的(波导)缝隙天线分析与设计
January,2008
大连海事大学学位论文原创性声明和使用授权说明
原创性声明
本人郑重声明：本论文是在导师的指导下,独立进行研究工作所取得的成果，撰写成硕士学位论文“基于时域有限差分法的缝隙天线分析与设计”。除论文中已经注明引用的内容外，中以明确方式标明.本论文中不包含任何未加明确注明的其他个人或集体已经公开发表或未公开发表的成果。
分类号
密级
U D C
单位代码10151
基于时域有限差分法的缝隙天线分析与设计
指导教师
房少军
职称
教授
学位授予单位
大连海事大学
申请学位级别
硕士
学科与专业
通信与信息系统
论文完成日期
论文答辩日期
答辩委员会主席
TheAnalysis and Design ofSlotAntenna Based onFinite-Difference Time-DomainMethod
Excited from a single end of a waveguide,the length of a slotted-waveguide antenna is restricted by the width of radar transmitted pulse, soit is difficult to implement alarge-scaleradarslot antenna with narrow beams. However, we should not ignore the fact that slot antenna has many advantages. To overcome this contradiction,the paper provides a design of slotted-waveguide antenna which is excited from both end of thewaveguide, and simulatesit with the help of HFSS. According to the requirement of directivity, this researcher designs the amplitudedistribution of the aperturesbased on modifiedChebyshev array, and calculatesit by writing program usingMATLAB. The complex arrangement of the slots in the narrow face of a waveguide brings much trouble to model it in HFSS, so the present paper offers a method to create the model byproducingmodelscriptsinMATLAB.This method can establish a complex modelquickly and exactly.The simulating result shows that this antenna can basically satisfy the practical requirement.
一种处理色散介质问题的通用时域有限差分方法
一种处理色散介质问题的通用时域有限差分方法3魏 兵 葛德彪 王 飞(西安电子科技大学物理系,西安 710071)(2007年12月17日收到;2008年4月11日收到修改稿) 色散介质的介电系数是频率的函数,使本构关系在时域成为卷积关系.这就给用时域有限差分方法计算色散介质中波的散射和传播带来了困难.现有算法往往要针对不同色散介质模型推导相应的递推公式,算法的通用性较差.本文完善和发展了移位算子2时域有限差分方法,使之成为一种处理色散介质电磁问题的通用方法.首先,证明了常见的三种色散介质模型(德拜模型、洛伦兹模型和德鲁模型)的介电系数均可以写成适于移位算子法计算的有理分式函数形式.然后,用9Π9t 代替j ω,过渡到时域,再引入时域移位算子z t 代替时间微分算子来处理有理分式函数形式的介电系数,给出离散时域本构关系的表示式,进而导出时域有限差分方法当中电位移矢量和电场强度之间的关系.最后,计算了几种色散介质的电磁散射,数值结果表明了本文方法和程序的通用性和正确有效性.关键词:时域有限差分方法,色散介质,移位算子PACC :4110H ,5170,52103国家自然科学基金(批准号:60871070)和国家博士后科学基金(批准号:20070421109)资助的课题. E 2mail :bwei @11引言近十几年来,随着国内外对时域有限差分(finite difference time domain ,FDT D )方法研究的深入,将该方法用于处理色散介质电磁问题引起人们的关注.1990年,Luebbers 等人[1]提出了适用于Debye 模型的递归卷积FDT D 方法(recursive conv olution FDT D ,RC 2FDT D ),然后将该方法推广到等离子体介质[2]和N 阶色散介质[3,4].Hunsberger 等人[5]将RC 2FDT D 方法推广用于磁化等离子体介质.Luebbers 等人[6]研究了色散介质球的电磁散射问题.K elley 等人[7]用电场的分段线性(piecewise linear )近似(P LRC 2FDT D 方法改善了RC 2FDT D 方法)的计算精度.Siushansian 等人[8]采用离散的梯形递归卷积(TRC 2FDT D )方法改善了RC 2FDT D 方法的计算精度.此外,处理色散介质电磁问题的FDT D 方法还有辅助方程(ADE )法[9—11],Z 变换法[12—14],电流密度卷积(J EC )法[15],Y oung 氏直接积分法[16—18],分段线性电流密度卷积(P LJ ERC )算法[19,20]等.近年来,对色散介质的研究已逐步深入到各向异性介质的情形[21—24].上述几种方法中,RC 法将电位移矢量写成电场强度的卷积,离散该卷积成迭代求和式,再联立电场强度和磁场强度的迭代式,实现时域迭代计算.J EC 法将极化电流密度表示为电场强度的卷积并离散得到迭代方程,再联立电场强度和磁场强度的迭代式,实现时域迭代计算.P LRC 法和P LJ ERC 法分别在RC 法和J EC 法的基础上引进分段线性近似以改善计算精度.ADE 法将Maxwell 方程和介质所满足的相关方程直接差分,得到一个包含多个量的差分方程组,从而实现场量的时域迭代计算.Z 变换法把频域本构关系变换到Z 域,然后再通过Z 域得到时域递推式.RC 法,P LRC 法,J EC 法和P LJ ERC 法等需进行复杂卷积计算.ADE 法和Z 变换法的数学过程也比较繁琐.总的来讲,现有方法往往需要对不同的色散介质推导相应的递推公式,并编制相应计算程序,算法和程序的通用性较差.2002年,葛德彪等人[25]提出了处理色散介质电磁问题的移位算子2时域有限差分(shift operator finite difference time domain ,S O 2FDT D )方法并讨论了该算法在非磁性等离子体中的应用.本文完善和发展了文献[25]所提出的S O 2FDT D 方法,使之成为色散介质电磁问题处理的通用算法.首先,证明了常见的三种色散介质模型(德拜模型、洛伦兹模型和德鲁模第57卷第10期2008年10月100023290Π2008Π57(10)Π6290208物 理 学 报ACT A PHY SIC A SI NIC AV ol.57,N o.10,October ,2008ν2008Chin.Phys.S oc.型)的介电系数均可以写成适于移位算子法计算的有理分式函数形式.然后,用9Π9t 代替j ω,过渡到时域,再引入时域移位算子代替时间微分算子来处理有理分式函数形式的介电系数.给出了用移位算子表达式代替对时间偏微分时从低阶到高阶的严格证明,推导了适应于三种色散介质模型的通用FDT D 递推表达式.最后,给出了本文方法的计算结果和文献结果以及解析解的比较.本文中时谐因子取exp (j ωt ).21常见色散介质模型及其介电系数表达式 常见的三种线性、各向同性色散介质模型的介电常数表达式[26]如下:1)德拜模型(Debye m odel )ε(ω)=ε∞+6pp =1εs ,p -ε∞,p1+j ωτp ≡ε∞+6pp =1Δεp1+j ωτp,(1)式中,Δεp =εs ,p -ε∞,p ,其中,εs ,p 为静态或零频时的相对介电系数,ε∞,p 为无穷大频率时的相对介电系数;τp 为极点弛豫时间.2)洛仑兹模型(Lorentz m odel )ε(ω)=ε∞+6pp =1Δεpω2p ω2p +2j ωδp -ω2,(2)式中Δεp 含义同(1)式,ωp 为极点频率,δp 为阻尼系数.3)德鲁模型(Drude m odel )ε(ω)=ε∞-6pp =1ω2pω2-j ωγp ,(3)式中ωp 为德鲁极点频率,γp 为极点弛豫时间的倒数.可以证明(见附录A ),上述几种色散介质模型中的相对介电系数εr (ω)均可以写成以下有理分式函数形式,即εr (ω)=6Nn =0p n (j ω)nΠ6Nn =0q n (j ω)n.(4)31各向同性色散介质FDT D 递推式 设各向同性色散介质频域中本构关系为D (ω)=ε(ω)E (ω),(5)式中ε(ω)为复数介电系数.色散介质中麦克斯韦旋度方程在无源情况为Δ×E =-9B Π9t =-μ9H Π9t ,Δ×H =9D Π9t .(6)(6)式中设导磁系数与频率无关,即B =μH .按照Y ee 元胞由(6)式第一式可得到磁场强度H 各分量随时间推进的计算公式[27].另外,将(6)式第二式在t =(n +1Π2)Δt 时刻作差分离散可得(Δ×H )n +1Π2=[D n +1-D n]ΠΔt ,(7)即Dn +1=D n +Δt (Δ×H )n +1Π2.(8)取电位移矢量D 的采样点与电场强度E 的采样点相同,离散(8)式可得D 各分量随时间推进计算公式为(以x 分量为例)D n +1x (i +1Π2,j ,k )=D nx (i +1Π2,j ,k )+Δt ・{[H n +1Π2z (i +1Π2,j +1Π2,k )-H n +1Π2z(i +1Π2,j -1Π2,k )]}ΠΔy -[H n +1Π2y(i +1Π2,j ,k +1Π2)-H n +1Π2y(i +1Π2,j ,k -1Π2)]ΠΔz }.(9)(9)式是从H 到D 的递推式,要使FDT D 计算进行下去,还需要知道从D 到E 的递推式.下面用移位算子法得到从D 到E 的递推式.41用移位算子法处理色散介质的本构关系 下面用移位算子法得到色散介质时域本构关系表达式,进而得到FDT D 所需的时域递推计算式.4111含移位算子的本构关系 设频域中介质的本构关系为(以x 分量为例)D x =ε0εr (ω)E x .(10)若介电系数εr (ω)可以写成以j ω为自变量的分式多项式,即(4)式的形式,利用频域到时域的算子转换关系j ω→9Π9t ,将(4)式代入(10)式得到时域本构关系为D x (t )=ε0εr (9Π9t )E x (t ),(11)式中εr (9Π9t )为相对介电系数的时域算子形式,表示为εr (9Π9t )=6Nn =0p n (9Π9t )nΠ6Nn =0q n (9Π9t )n.(12)192610期魏 兵等:一种处理色散介质问题的通用时域有限差分方法将(12)式代入(11)式,并将分母上的求导运算移到等式左边可得6Nn =0q n (9Π9t )nD x (t )=ε06Nn =0p n (9Π9t )nE x (t ).(13)(13)式是时域中含时间导数算子的本构关系,这是一个微分方程.为得到方程(13)在时域的递推计算式,下面讨论时间导数算子在离散时域中的形式.设函数y (t )=9f (t )Π9t ,(14)(14)式在(n +1Π2)Δt 的中心差分近似为(y n +1+y n )Π2=(f n +1-f n)ΠΔt ,(15)其中(15)式左端取平均值近似.引进离散时域的移位算子z t ,定义为z t fn=fn +1,(16)即移位算子的作用相当于使离散时域数列的n 时刻值移位到函数在n +1时刻的值.合并(15)和(16)式可得y n=(2Δt )[(z t -1)Π(z t +1)]f n ,(17)比较(14)和(17)式有9Π9t →(2ΠΔt )[(z t -1)Π(z t +1)].(18)可以证明,高阶时间导数的移位算子形式为(见附录B )(9Π9t )n→{(2ΠΔt )[(z t -1)Π(z t +1)]}n .(19)将(19)式代入(13)式并整理得到以下离散时域的本构关系(为简化,以下令h =2ΠΔt ):6Nl =0q l {h [(z t -1)Π(z t +1)]}lD nx=ε06Nl =0p l {h [(z t -1)Π(z t +1)]}lE nx .(20)将(20)式两边乘以(z t +1)N得6Nl =0q l h l (z t +1)N -l(z t -1)l D nx=ε06Nl =0p l h l (z t +1)N -l(z t -1)l E n x .(21)(21)式称为离散时域含移位算子的本构关系.4121离散时域的递推关系 首先考虑(2)式中有理分式表示中的最高幂次N =1和N =2的情形.情形1 令N =1,将N 代入(21)式整理得[(q 0+q 1h )z t +(q 0-q 1h )]D nx=ε0[(p 0+p 1h )z t +(p 0-p 1h )]E nx .(22)根据(16)式,可将(22)式写为E n +1x =[a 0(D n +1x Πε0)+a 1(D nx Πε0)-b 1E nx ]Πb 0.(23)(23)式给出了从D 到E 的递推计算公式,其中a 0=q 0+q 1h ,a 1=q 0-q 1h ,b 0=p 0+p 1h ,b 1=p 0-p 1h .(24) 情形2 令N =2,将N 代入(21)式整理得{[q 0+q 1h +q 2h 2]z 2t +[2q 0-2q 2h 2]z t +[q 0-q 1h +q 2h 2]}D nx={[p 0+p 1h +p 2h 2]z 2t +[2p 0-2p 2h 2]z t+[p 0-p 1h +p 2h 2]}ε0E nx .(25)根据(16)和(21)式,(25)式可写为E n +1x =[a 0(D n +1x Πε0)+a 1(D nx Πε0)+a 2(D n -1x Πε0)-b 1E n x -b 2E n -1x]Πb 0. (26)(26)式给出了从D 到E 的递推计算公式,其中a 0=q 0+q 1h +q 2h 2,a 1=2q 0-2q 2h 2,a 2=q 0-q 1h +q 2h 2,b 0=p 0+p 1h +p 2h 2,b 1=2p 0-2p 2h 2,b 2=p 0-p 1h +p 2h 2.(27) 对于N ≥3的一般情况,(21)式可写为En +1x=1b 06Nl =0a l (Dn +1-l xΠε0)-6Nl =1b l E n +1-lx.(28)式中a l ,b l 可由q 0,q 1,…,q N ,p 0,p 1,…,p N 表示.可以看出,当色散介质的相对介电系数的分式有理函数的最高次幂N =1时,只需知道前一时刻的D 与E 值和当前时刻的D 值即可求出当前时刻的E .而当N =2时,则需要知道前两个时刻的D ,E 值和当前时刻的D 值才能求出当前时刻的E .51数值结果 算例中入射波均取高斯脉冲E i (t )=exp [-4π(t -t 0)2Πτ2],(29)式中τ为脉冲宽度.以下δ均表示FDT D 计算时的元胞尺度,Δt 表示时间步长,取Δt =δΠ(2c ),其中c 为真空光速.算例1 吸波介质球的散射吸波介质球半径为0125米,非磁性介质(μr =2926物 理 学 报57卷110,σm =0).其复相对介电常数可以表示为ε(ω)=ε∞+(εs -ε∞)Π(1+j ωt 0)-j σΠωε0,(30)式中,ω为入射波频率,j 为虚数单位,εs =1116为频率为零时的相对介电常数,ε∞=1101为频率为无穷大时的相对介电常数,电导率σ=2195×10-4Ω,t 0=41497×10-10s.图1是该介质球的后向雷达散射截面(RCS ),图中实线是考虑色散现象后用色散介质FDT D 程序的计算结果.圆点是Mie 级数的计算结果.可以看出,色散介质FDT D 的计算结果与Mie级数的结果符合的很好.计算中δ=313×10-3m ,τ=60Δt.图1 吸波球的后向RCS算例2 非磁性等离子体球的后向RCS 非磁化等离子体的相对介电系数[28]为εr (ω)=[1+ω2p Π(ω(j v c -ω))]={[(j ω)2+v c (j ω)+ω2p ]Π[(j ω)2+v c (j ω)]},(31)式中,ω2p =(Ne 2)Π(m ε0)为等离子体频率(N ,e ,m和ε0分别为电子密度、电子电量、电子质量和真空介电常数),v c 为电子平均碰撞频率.等离子体球半径为3175mm ,非磁性等离子体的相对介电常数可以表示为(31)式.取v c =210×1010s -1,ωp =118×1011H z ,δ=510×10-2mm ,τ=60Δt .图2和图3分别为等离子体后向散射场的时域波形和后向RCS.图中圆圈表示由Mie 级数所得的结果,实线表示本文方法的计算结果,可以看出两者相符.算例3 等离子体覆盖导体圆柱的宽频后向散射特性覆盖等离子体的导体圆柱的截面如图4所示.图2 等离子体球后向散射场时域波形图3 等离子体球的后向RCS图4 等离子体覆盖导体圆柱的示意图导体圆柱半径为a =012m ,d =b -a 表示等离子体包层厚度.入射方向如图,后向散射.计算中取d =011m ,δ=015mm.T M 波入射时柱体的散射宽度如图5所示.图5(a )中实线为本文方法的计算结果,○和★是文献[29]的结果,两者符合很好.等离子体的浓度为N =510×1017m -3.由图可见,等离子体对392610期魏 兵等:一种处理色散介质问题的通用时域有限差分方法图5 二维覆盖等离子体的金属圆柱的散射宽度 (a )T M 情形;(b )TE 情形电磁波的吸收效果明显,且具有很宽的吸收频带.61结论本文的推证和计算结果说明,S O 2FDT D 是一种普遍适用于三种常见色散介质模型的通用计算方法.该方法数学推证过程清晰明了,便于编制通用程序处理不同类型色散介质的电磁散射问题,且具有较小的内存需求量.本文的数值结果说明了移位算子FDT D 方法的正确性和通用性.附录A 三种色散介质模型的有理分式函数形式 S O 2FDT D 方法要求,色散介质的相对介电系数可以写成以j ω为自变量的有理分式函数形式.由于介电系数和相对介电系数只差一个常量ε0,下面用归纳法证明三种常见色散介质模型的介电系数可以写成以j ω为自变量的有理分式函数形式,即ε(ω)=6Nn =0p n (j ω)nΠ6Nn =0q n (j ω)n.(A1) 1)德拜模型(Debye model)ε(ω)=ε∞+6pp =1[Δεp Π(1+j ωτp )].(A2) 证明 当P =1时ε1(ω)=[(ε∞+Δε1(j ω)0+τ1ε∞(j ω)1]Π[(j ω)0+τ1(j ω)1],(A3)可见,当p =1时,(A2)式可以写成(A1)式的形式.设当p =n (n ≥1)时,(A2)式可以写成(A1)式的形式,即εn (ω)=ε∞+6np =1Δεn 1+j ωτp=ε06mn =0p n (jω)nΠ6mn =0q n (j ω)n,(A4)则当p =n +1时,εn +1(ω)=(p 0b +q 0a )(j ω)0+(p 1b +p 0c +q 1a )(j ω)1q 0b (j ω)0+(q 1b +q 0c )(j ω)1+(q 2b +q 1c )(j ω)2+…×+(p 2b +p 1c +q 2a )(j ω)2+…+(p m b +p m -1c +q m a )(j ω)m +p m c (j ω)m +1+(q m b +q m -1c )(j ω)m +q m c (j ω)m +1,(A5)上式中a =Δεn +1(j ω)0Πε0,b =(j ω)0,c =τn +1.故εn +1(ω)也可以写成6m +1n =0p n (j ω)nΠ6m +1n =0q n (j ω)n的形式.综上所述,无论p 取何值,(A2)式都可以写成(A1)式的形式. 2)洛仑兹模型(Lorentz model)ε(ω)=ε∞+6pp =1Δεpω2pΠ(ω2p +2j ωδp -ω2).(A6)4926物 理 学 报57卷 证明 当p =1时ε1(ω)=(ε∞ω21+Δε1ω21)(j ω)0+2ε∞δ(j ω1)1+ε∞(j ω)2ω21(j ω)0+2δ(j ω)1+(j ω)2,(A7)故,当p =1时,(A6)式以写成(A1)式的形式.设当p =n (n ≥1)时,(A6)式可以写成(A1)式的形式,即εn (ω)=ε∞+6np =1Δεpω2p ω2p +2j ωδp -ω2=6mn =0p n (j ω)nΠ6mn =0q n (j ω)n,(A8)则当p =n +1时εn +1(ω)=(p 0b +q 0a )(j ω)0+(p 1b +p 0c +q 1a )(j ω)1+(p 2b +p 1c +p 0d +q 2a )(j ω)2…q 0b (j ω)0+(q 1b +q 0c )(j ω)1+(q 2b +q 1c +q 0d )(j ω)2+…×+(p m b +p m -1c +p m -2d +q m a )(j ω)m +(p m c +p m -1d )(j ω)m +1+p m d (j ω)m +2+(q m b +q m -1c +q m -2d )(j ω)m +(q m c +q m -1d )(j ω)m +1+q m d (j ω)m +2,(A9)上式中a =Δεn +1ω2n +1Πε0,b =ω2n +1,c =2δn +1,d =1,故εn +1(ω)也可以写成6m +2n =0p n (j ω)nΠ6m +2n =0q n (j ω)n的形式.综上所述,无论p 取何值,(A6)式都可以写成(A1)式的形式.3)德鲁模型(Drude model)ε(ω)=ε∞-6pp =1ω2pΠ(ω2p -j ωγp ).(A10) 证明 当p =1时ε1(ω)=[(ε∞ω21-ω21)(j ω)0-ε∞γ1(j ω)1]Π[ω21(j ω)0-γ1(j ω)1],(A11)显然,当p =1时,(A10)式可以写成(A1)式的形式.设当p =n(n ≥1)时,(A10)式可以写成(A1)式的形式,即εn (ω)=6mn =0p n (j ω)nΠ6mn =0q n (j ω)n,(A12)则当p =n +1时εn +1(ω)=(p 0b +q 0a )(j ω)0+(p 1b +p 0c +q 1a )(j ω)1q 0b (j ω)0+(q 1b +q 0c )(j ω)1+(q 2b +q 1c )(j ω)2+…×+(p 2b +p 1c +q 2a )(j ω)2+…+(p m b +p m -1c +q m a )(j ω)m +p m c (j ω)m +1+(q m b +q m -1c )(j ω)m +q m c (j ω)m +1,(A13)上式中a =-ω2n +1(j ω)0Πε0,b =ω2n +1(j ω)0,c =-γn +1,所以εn +1(ω)也可以写成6m +1n =0p n (j ω)nΠ6m +1n =0q n (jω)n的形式.综上所述,无论p 取何值,(A10)式均可以写成(A1)式的形式.附录B 高阶时间导数移位算子形式的证明已知9Π9t →(2ΠΔt )[(z t -1)Π(z t +1)],以下用归纳法证明(9Π9t )n→{(2ΠΔt )[(z t -1)Π(z t +1)]}n . 证明 设函数y 为函数f 对时间的二阶导数y (t )=92f (t )Π9t 2,(B1)将上式右端在(n +1Π2)Δt 的中心差分近似为yn=(fn +1-2fn+fn -1)ΠΔt 2,(B2)将(B2)式左端函数y 在n 时刻的值由n +1时刻,n 时刻和n-1时刻的值平均得到,则(B2)式可以改写为yn=y n2+yn2=yn +1+yn -14+yn2=yn +1+2y n+y n -14=fn +1-2f n+f n -1Δt 2.(B 3) 引进离散时域的移位算子z t ,即定义为(16)式,由(16)和(B3)式可得592610期魏 兵等:一种处理色散介质问题的通用时域有限差分方法y n-1=4Δt2z2t-2z t+1z2t+2z t+1f n-1=2Δtz t-1z t+12f n-1.(B4)可见当n=2时(19)式成立.设n=k时(19)式成立,即(9Π9t)k→{(2ΠΔt)[(z t-1)Π(z t+1)]}k,(B5)则当n=k+1时,设函数y是函数f对时间的k+1阶偏导y(t)=9k+1f(t)Π9t k+1=(9Π9t)[9k f(t)]Π9t k,(B6)令x(t)=9k f(t)Π9t k,(B7) (B6)式可以写成y(t)=(9Π9t)(x(t)).(B8)同样引进离散时域的移位算子zt,即z t x n=x n+1,(B9)可得y n=(2ΠΔt)[(z t-1)Π(z t+1)]x n.(B10)综合(B5),(B7)和(B10)式可得y k+1=2Δtz t-1z t+12Δtz t-1z t+1k=2Δtz t-1z t+1k+1.(B11)可见当n=k+1时(19)式是成立的.综上所述,无论k取任何整数值,(19)式均成立.[1]Luebbers R J,Hunsberger F,K unz K S,S tandler R,Schneier M1990IEEE pat.32222[2]Luebbers R J,Hunsberger F,K unz K S1991IEEE Trans.AntennasPropagat.3929[3]Luebbers R J,Hunsberger F R1992IEEE Trans.AntennasPropagat.401297[4]P ontalti R,Cristoforetti L,Antolini R,Cescatti L1994IEEE Trans.Microwave Theory Tech.42526[5]Hunsberger F,Lubbers R J,K unz K S1992IEEE Trans.AntennasPropagat.401489[6]Luebbers R J,S teich D,K unz K1993IEEE Trans.AntennasPropagat.411249[7]K elley D F,Luebbers R J1996IEEE Trans.Antennas Propagat.44792[8]S iushansian R,Lovetri J1995IEEE Microwave Guided Wave Lett.5426[9]Nickisch L J,Franke P M1992IEEE Antennas Propagat.Mag.3433[10]G andhi O P,G ao B Q,Chen T Y1993IEEE Trans.MicrowaveTech.41658[11]T akayama Y,K laus W1994IEEE Microw.Wireless Compon.Lett.12102[12]Sullivan D M1992IEEE Trans.Antennas Propagat.401223[13]Sullivan D M1995IEEE Trans.Antennas Propagat.43676[14]Sullivan D M1996IEEE Trans.Antennas Propagat.4428[15]Chen Q,K atsurai M,A oyagi P H1998IEEE Trans.AntennasPropagat.461739[16]Y oung J L1994Radio Sci.291513[17]Y oung J L1995IEEE Trans.Antennas Propagat.43422[18]Y oung J L1996IEEE Trans.Antennas Propagat.441283[19]Xu L J,Y uan N C2005IEEE microwave and Wireless ComponentsLetter s15277[20]Liu S B,M o J J,Y uan N C2004Acta Phys.Sin.53778(inChinese)[刘少斌、莫锦军、袁乃昌2004物理学报53778] [21]Y ang L X,G e D B,W ei B2007Acta Phys.Sin.564509(inChinese)[杨利霞、葛德彪、魏 兵2007物理学报564509] [22]Liu S B,M o J J,Y uan N C2004Acta Phys.Sin.53783(inChinese)[刘少斌、莫锦军、袁乃昌2004物理学报53783] [23]Liu S B,M o J J,Y uan N C2004Acta Phys.Sin.532233(inChinese)[刘少斌、莫锦军、袁乃昌2004物理学报532233] [24]Y ang L X,G e D B2006Acta Phys.Sin.551751(in Chinese)[杨利霞、葛德彪2006物理学报551751][25]G e D B,Wu YL,Zhu X Q2003Chin.J.Radio Sci.18359(inChinese)[葛德彪、吴跃丽、朱湘琴2003电波科学学报18359][26]Allen T,Susan C H2005Computational Electrodynamics2the FiniteDifference Time Domain Method(third edition)(London:ArtechH ouse)[27]G e D B,Y an Y B2005Finite2Difference Time2Domain Method forElectromagnetic Waves(2thed)(X iπan:X idian University Press)(inChinese)[葛德彪、闫玉波2005电磁波时域有限差分方法(第二版)(西安:西安电子科技大学出版社)][28]G inzburg V L1970The Propagation o f Electromagnetic Waves inPlasmas(2thed)(Ox ford:Pergam on Press)[29]M o J J,Liu S B,Y uan N C2003J.Microwaves1920(in Chinese)[莫锦军、刘少斌、袁乃昌2004微波学报1920]6926物 理 学 报57卷A general method for finite difference time domain modelingof wave propagation in frequency 2dispersive media 3W ei Bing G e De 2Biao W ang Fei(Department o f Physics ,Xi πdian Univer sity ,Xi πan 710071,China )(Received 17December 2007;revised manuscript received 11April 2008)AbstractThe analysis of electromagnetic scattering and propagation in dispersive media is com plicated in time domain ,because its dielectric property is frequency 2dependent.A disadvantage of the prevailing alg orithms is the need to deduce different formulations for each dispersion m odel.In this paper ,the shift operator finite difference time domain (S O 2FDT D )method is developed.First ,we prove that the com plex perm ittivity of three kinds of general dispersive media m odels ,i.e.Debye m odel ,the Lorentz m odel and the Drude m odel ,may be described by rational polynom ial functions in j ω.By introducing a shift operator z t ,the constitutive relation between D and E is derived in discretised time domain.The shift operator method is then applied to the general dispersive medium case.The recursive formulation for D and E available for FDT D com putation is obtained.Finally ,the scatterings by a dispersive sphere and a PEC object covered w ith dispersive media are com puted.The com puted results are in g ood agreement w ith the literature and the one obtained by M ie πs series solution.This illustrates the generality and the feasibility of the presented scheme.K eyw ords :FDT D method ,dispersive media ,shift operator PACC :4110H ,5170,52103Project supported by the National Natural Science F oundation of China (G rant N o.60871070)and the National Science F oundation for P ost 2doctoralScientists of China (G rant N o.20070421109).E 2mail :bwei @792610期魏 兵等:一种处理色散介质问题的通用时域有限差分方法。
基于时域有限差分法的gis设备暂态电场仿真计算
摘要由于气体绝缘变电站(Gas Insulated Substation)设备结构紧凑，遭受过电压时其幅值大、频率高，而SF6的绝缘特性不光取决于绝缘结构中电场的不均匀程度，还与作用电压源的陡度密切相关，因此，研究GIS设备暂态电场分布对指导绝缘设计十分必要。
由于时域有限差分法(FDTD)特有的优越性，本文采用FDTD方法对GIS中的暂态电场分布进行仿真研究。
首先，为提高将FDTD计算GIS电场时算法的性能，对于算法的各个环节进行了探讨，得出研究的重点在于选择合适的网格剖分算法和吸收边界条件；随后，按照建模时的结构特点将GIS分为同轴母线管道和GIS复杂设备两部分，对其中的网格剖分算法等关键性问题进行研究。
为实现对GIS同轴母线管道的建模，提出了一套导体-介质体通用的共形网格技术，并对共形网格信息生成等关键技术进行了研究，然后对共形网格技术的稳定性进行分析，比较了不同因素对稳定性的影响。
为实现复杂设备的准确计算并提高计算效率，采用了非均匀网格方法。
对复杂设备剖分时的新问题进行了研究，解决了介质存在重叠和缝隙时共形网格信息的生成问题，提出了一种新的奇异点的处理方法。
对算法有效性进行了验证，并给出了GIS设备计算时吸收边界条件的选取规则。
采用上述FDTD方法，以上海思源电气公司制造的500kV的GIS母线筒和隔离开关室为例，研究了其在特快速暂态过电压(VFTO)和雷电波下的暂态电场特性，验证了研究暂态电场的意义。
关键词：时域有限差分法GIS设备暂态电场计算共形网格技术AbstractThe over-voltage in GIS equipment usually has high amplitude and frequency due to the compact structure. The insulation characteristics of SF6gas not only depends on the inhomogeneity of electric field in insulation structure, but also are closely related to the steepness of the voltage source. Therefore, the research of the dynamic electric field distribution in GIS is very necessary to guide insulation design. Because of the superiority of time domain finite difference method (FDTD), the FDTD method is used to simulate the dynamic electric field distribution in GIS.First of all, the key problems in each procedure of the FDTD simulation are discussed to improve the performance in the GIS electric field calculation. We find that the focal points include two aspects, i.e., gridding algorithm and absorbing boundary condition. Then, we classify the GIS into coaxial bus duct and GIS equipment according to its structure characteristic, and study the key problems.In order to model the GIS coaxial bus duct, a general technology which is used to gridding the conductor and dielectric is put forward, and the key technologies of conformal grid information generation are introduced. Then we analysis the stability of the gridding technology and compare the influence of different factors. In order to calculate the electric field of complex equipment accurately and improve the computational efficiency, the non-uniform gridding method is introduced. New problems in modeling complex equipment are studied, including the generation of conformal mesh information and how to deal with singular points. Then the validity of the algorithm is verified, and the selection rules of the absorption boundary condition are given.Taking the 500kV GIS bus duct and disconnector made by Shanghai Siyuan electric company as study example, the transient electric field under very fast transient overvoltage (VFTO) and lightning wave is calculated using FDTD method, which validates the significance of calculating transient electric field.Keywords: Time domain finite difference method (FDTD) GIS equipmentTransient electric field calculation Conformal grid technology目录摘要 (I)Abstract .............................................................................................................. I I 1绪论 (1)1.1研究背景和意义 (1)1.2国内外研究现状 (3)1.3本文的主要工作 (8)2GIS同轴母线管道建模方法研究 (10)2.1FDTD共形技术 (10)2.2共形网格自动剖分算法 (20)2.3算法有效性验证 (26)2.4本章小结 (30)3复杂设备的建模方法 (31)3.1非均匀网格算法 (31)3.2复杂设备中的剖分新问题 (36)3.3算法有效性验证 (44)3.4本章小结 (48)4GIS设备中的暂态电场分布 (49)4.1不同电压下GIS母线管道中的暂态电场分布 (49)4.2不同电压下隔离开关室中的暂态电场分布 (59)4.3本章小结 (65)5总结和展望 (66)5.1全文总结 (66)5.2展望 (68)致谢 (69)参考文献 (70)附录1 攻读学位期间发表的论文及参与的科研项目 (78)1绪论1.1研究背景和意义1.1.1GIS简介GIS是气体绝缘变电站(Gas Insulated Substation)的简称[1]，它将断路器、隔离开关、电压互感器、电流互感器、母线、进出线套管等设备封闭在金属壳内，主要采用SF6气体作为绝缘和灭弧介质。
时域有限差分法
Problem 5.1In this illustrative solution, the electric-field hard source condition of (5.1) is implemented at the far-left grid boundary. The source time function has an amplitude of 1.0 V/m and a frequency of 10 GHz. The reflecting barrier (PEC) is implemented at the far-right grid boundary. The computational domain represents a physical length of 15 cm.Matlab code:%***********************************************************************% 1D FINITE-DIFFERENCE TIME-DOMAIN SOLUTION: PLANE WAVE PROPAGATION%***********************************************************************%% Program author:% Prof. Susan C. Hagness% Department of Electrical and Computer Engineering% University of Wisconsin-Madison% 1415 Engineering Drive% Madison, WI 53706-1691% hagness@%%***********************************************************************clear;%..........Material Parameters............cc=2.99792458e8; %speed of light in free spacemuz=4.0*pi*1.0e-7; %permeability of free spaceepsz=1.0/(cc*cc*muz); %permittivity of free spaceeps=[1.0]; %relative permittivitysig=[0.0]; %electric conductivitymur=[1.0]; %relative permeabilitysim=[0.0]; %magnetic lossmedia=length(eps);%..........Space, Time, and Source Parameters...S=1.0;freq=10e9; %frequency of sinusoidal excitation = 10 GHzE0=1.0; %amplitude of sinusoidal excitation = 1.0 V/mlambda=cc/freq;length=0.15; %physical length of grid (in units of m)dx=lambda/20; %grid resolution of 20 cells per wavelengthdt=S*dx/cc;ie=round(length/dx)+1; %number of Ez samples in gridih=ie-1; %number of Hy samples in gridnmax=3*round(ie*S);source(1:nmax)=E0*sin(2*pi*freq*(1:nmax)*dt);%..........Initial Conditions...........ez(1:ie)=0.0;hy(1:ih)=0.0;Problem 5.1 (cont.)%..........Update Coefficients.........for i=1:mediaeprop=sig(i)*dt/(2.0*epsz*eps(i));ca(i)=(1.0-eprop)/(1.0+eprop);cb(i)=dt/(epsz*eps(i)*dx)/(1.0+eprop);hprop=sim(i)*dt/(2.0*muz*mur(i));da(i)=(1.0-hprop)/(1.0+hprop);db(i)=dt/(muz*mur(i)*dx)/(1.0+hprop);endmediaez(1:ie)=1; %media pointersmediahy(1:ih)=1;%..........Time-Stepping..............for n=1:nmax%%%%%%%%%%% EZ UPDATES %%%%%%%%%%%ez(2:ih)=ca(mediaez(2:ih)).*ez(2:ih)+cb(mediaez(2:ih)).*(hy(2:ih)-hy(1:ih-1)); ez(1)=source(n); % incident wave source condition (hard source) %%%%%%%%%%% HY UPDATES %%%%%%%%%%%hy(1:ih)=da(mediahy(1:ih)).*hy(1:ih)+db(mediahy(1:ih)).*(ez(2:ie)-ez(1:ih));%%%%%%%%%%% DATA OUTPUT %%%%%%%%%%%subplot(2,1,1),plot(ez,'r'),axis([1 ie -3 3]);ylabel('Ez'); title(['n=',num2str(n)]);subplot(2,1,2),plot(hy,'b'),axis([1 ih -8e-3 8e-3]);ylabel('Hy'); xlabel('grid coordinate')pause(0.05)endProblem 5.1 (cont.)Visualizations of the electric and magnetic field distributions:The incident sinusoidal wave has been launched from the hard source at the far-left grid boundary and is propagating in the +x direction. The amplitude of the traveling wave is 1.0 V/m.At time step n=100, the leading edge of theincident wave has reached the PEC mirror at the far-right grid boundary.The wave reflected from the PEC boundary propagates in the –x direction. A standing wave is established in the region of the grid where the reflected and incident waves are overlapping in space. The amplitude of the standing wave is 2.0 V/m. The remainder of the grid contains only the incident traveling wave of amplitude 1.0 V/m.At time step n=145, the leading edge of the reflected wave has reached grid coordinate i=55 (i.e. the wave has propagated a total distance of 100+45 grid cells in 145 time steps).Traveling waveStanding waveProblem 5.1 (cont.)At time step n=200, the reflected wave hasreached the far-left grid boundary. The entiregrid is now filled with the standing wave.The wave is retro-reflected from the hard sourceat the far-left grid boundary. A new standingwave is established in the region of the gridwhere the retro-reflected, reflected, and incidentwaves are counter-propagating. The amplitude ofthis new standing wave is 3.0 V/m.Problem 5.4The discrete-time parameters of the Gaussian pulse are chosen to replicate the example shown in Fig. 5.9. Since no physical dimensions are given in the context of Fig. 5.9, the choice of the grid cell size is an arbitrary one.Matlab code:%***********************************************************************% 1D FINITE-DIFFERENCE TIME-DOMAIN FORMULATION WITH% TOTAL-FIELD/SCATTERED-FIELD REGIONS%***********************************************************************%% Program author:% Prof. Susan C. Hagness% Department of Electrical and Computer Engineering% University of Wisconsin-Madison% 1415 Engineering Drive% Madison, WI 53706-1691% hagness@%%***********************************************************************clear;%..........Material Parameters............cc=2.99792458e8; %speed of light in free spacemuz=4.0*pi*1.0e-7; %permeability of free spaceepsz=1.0/(cc*cc*muz); %permittivity of free spaceeps=[1.0]; %relative permittivitysig=[0.0]; %electric conductivitymur=[1.0]; %relative permeabilitysim=[0.0]; %magnetic lossmedia=length(eps);%..........Space and Time Parameters...S=1.0;dx=1.0;dt=S*dx/cc;ie=400; %number of Ez samples in gridih=ie-1; %number of Hy samples in gridileft=100; %location of the left TF/SF boundaryiright=300; %location of the right TF/SF boundarynmax=500;%..........Source Parameters...........ntau=20;n0=3*ntau;source(1:nmax)=exp(-((1:nmax)-n0).^2/ntau^2);Problem 5.4 (cont.)%..........Initial conditions..........ez(1:ie)=0.0;hy(1:ih)=0.0;ezi(1:ie)=0.0; %incident-wave gridhyi(1:ih)=0.0; %incident-wave grid%..........Update Coefficients.........for i=1:mediaeprop=sig(i)*dt/(2.0*epsz*eps(i));ca(i)=(1.0-eprop)/(1.0+eprop);cb(i)=dt/(epsz*eps(i)*dx)/(1.0+eprop);hprop=sim(i)*dt/(2.0*muz*mur(i));da(i)=(1.0-hprop)/(1.0+hprop);db(i)=dt/(muz*mur(i)*dx)/(1.0+hprop);endcaez(1:ie)=ca(1);cbez(1:ie)=cb(1);dahy(1:ih)=da(1);dbhy(1:ih)=db(1);PEC=input('PEC mirror? ');if PEC=='y'caez(ie/2)=-1.0;cbez(ie/2)=0.0;end%..........Time-Stepping..............for n=1:nmax%%%%%%%%%%% EZINC UPDATES %%%%%%%%%%%ezi(1)=source(n);rbcir=ezi(ih);ezi(2:ih)=ca(1)*ezi(2:ih)+cb(1)*(hyi(2:ih)-hyi(1:ih-1));ezi(ie)=rbcir; %simple radiation boundary condition (right boundary) %%%%%%%%%%% EZ UPDATES %%%%%%%%%%%rbcl=ez(2);rbcr=ez(ih);ez(2:ih)=caez(2:ih).*ez(2:ih)+cbez(2:ih).*(hy(2:ih)-hy(1:ih-1));ez(1)=rbcl; %simple radiation boundary condition (left boundary) ez(ie)=rbcr; %simple radiation boundary condition (right boundary) % ...correction at TFSF interface...ez(ileft)=ez(ileft)-cbez(ileft)*hyi(ileft-1); %Eq. 5.28ez(iright)=ez(iright)+cbez(iright)*hyi(iright); %Eq. 5.34 %%%%%%%%%%% HYINC UPDATES %%%%%%%%%%%hyi(1:ih)=da(1)*hyi(1:ih) + db(1)*(ezi(2:ie)-ezi(1:ih));%%%%%%%%%%% HY UPDATES %%%%%%%%%%%hy(1:ih)=dahy(1:ih).*hy(1:ih)+dbhy(1:ih).*(ez(2:ie)-ez(1:ih));% ...correction at TFSF interface...hy(ileft-1)=hy(ileft-1)-dbhy(ileft-1)*ezi(ileft); %Eq. 5.32hy(iright)=hy(iright)+dbhy(iright)*ezi(iright); %Eq. 5.35Problem 5.4 (cont.)%%%%%%%%%%% DATA OUTPUT %%%%%%%%%%%plot(ez,'r-'),axis([1 ie -1 1]);xlabel('Grid i coordinate'); ylabel('Electric field E_z'); title(['n=',num2str(n)]);line([ileft ileft],[-1 1]); line([iright iright],[-1 1]);if PEC=='y'line([ie/2 ie/2],[0 1])endpause(0.05)endProblem 5.4 (cont.)Fig. 5.9(a) – free space everywhere in the gridProblem 5.4 (cont.)Fig. 5.9(b) – PEC mirror placed within the total-field region。
一种提高内存使用效率的时域有限差分算法
一种提高内存使用效率的时域有限差分算法张品;陈亦望;傅强【摘要】证明了即使在无源区域，局部一维时域有限差分法（LOD～FDTD）所给出的电磁场量也不满足零散度关系，推导了该散度关系的具体表达式。
基于该非零散度关系和麦克斯韦旋度方程，将LOD-FDTD法与减缩时域有限差分法（R—FDTD）相结合，得到一种新的局部一维减缩时域有限差分法（LOD-R—FDTD）。
该方法不仅具有LOD—FDTD方法的优势，计算公式简单，消除了CFL稳定条件对时间步长的限制，而且与LOD—FDTD相比平均节约了1／3内存使用量。
通过仿真计算与其他方法对比，证明了LOD～R—FDTD方法的准确性和有效性。
%In this paper, it is proven that the divergence relationship of electric-field and magnetic-field is non-zero even in charge-free regions when the electric-field and magnetic-field are calculated with locally one-dimensional finite-difference time-domain（LOD-FDTD） method, and the concrete expression of the divergence relationship is derived. Based on the non-zero divergence relationship and Maxwell curl e- quations set, the LOD-FDTD method which is unconditionally stable is combined with the reduced finite-difference time-dornain（R-FDTD）method. In the LOD-R- FDTD method, the advantage of LOD-FDTD is combined, which having simpler formulation to eliminate the restraint of the courant-friedrich-levy （CFL）condition. The memory requirement of LOD-R FDTD is reduced by 1/3 compared with LODFDTD in average. The formulation is presented and the accuracy and efficiency of the proposed method is verified by comparing the results with the conventional results.【期刊名称】《电波科学学报》【年(卷),期】2011(026)004【总页数】6页(P814-819)【关键词】时域有限差分方法;局部一维时域有限差分法;减缩时域有限差分法【作者】张品;陈亦望;傅强【作者单位】解放军理工大学工程兵工程学院,江苏南京210007;解放军理工大学工程兵工程学院,江苏南京210007;解放军理工大学工程兵工程学院,江苏南京210007【正文语种】中文【中图分类】O4411.引言时域有限差分方法（FDTD）已经广泛应用于电磁散射计算领域［1－3］。
时域有限差分方法
时域有限差分方法
《时域有限差分方法》
嘿，你知道吗，有一种超厉害的方法叫时域有限差分方法！这可真是个神奇的玩意儿。
想象一下，我们要研究那些看不见摸不着的电磁波啊之类的东西。
以前可麻烦了，但有了时域有限差分方法，就好像打开了一扇新的大门。
它是怎么工作的呢？简单来说，就是把我们要研究的区域划分成很多很多小格子，就像一个大拼图一样。
然后呢，通过计算这些小格子之间的变化，来了解整个区域的情况。
这个方法的好处可多啦！它能处理各种复杂的情况，不管是奇形怪状的物体，还是变化多端的环境。
而且，它很直观，让我们能清楚地看到电磁波是怎么传播、怎么变化的。
在实际应用中，时域有限差分方法可太有用了。
比如在通信领域，它能帮助我们设计更好的天线，让信号传输得更远更稳定。
在雷达系统中，它能让我们更准确地探测目标。
我觉得时域有限差分方法真的是一项非常了不起的技术，给我们探索和理解各种物理现象带来了巨大的帮助。
时域有限差分法及其在开关电源电磁辐射计算中的应用
时域有限差分法与其它电磁场数值计算方法相比具有许多突出的线路板（ＰＣＢ）的设计方面，与信号处理电路较规则的布线不同， 开关电源的线路板布线具有较大的随意性，由此增加了ＥＭｌ分析中ＰＣＢ的分布参
数提取的难度； ５、与信号处理电路中线路阻抗匹配的情况不同，开关电源的干扰源阻抗与网
络阻抗不仅不匹配，而且随工况变化，这无疑给ＥＭＩ滤波器的设计带来了一定的 困难。同时ＥＭＩ滤波器中的Ｌ、Ｃ元件还必须承受很大的无功功率，因而降低了 开关电源的整体效率，增大了开关电源的体积。
Ｉｎ ｔｈｅ ｓｔｕｄｙ Ｏｉｌ ｔｈｅ ｐｅｒｆｏｒｍａｎｃｅｓ ｏｆ ＰＭＬ ａｂｓｏｒｂｉｎｇ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎｓ，ｔｈｅ ｔｈｅｓｉｓ ｉｍｐｒｏｖｅｄ，ｒｅｓｐｅｃｔｉｖｅｌｙ，ｔｈｅ ２Ｄ ａｎｄ ３Ｄ ＦＤＴＤ ｍｅｔｈｏｄｓ ｗｉｔｈ ＰＭＬ ａｂｓｏｒｂｉｎｇ
ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ，ａｎｄ ｃｏｎｓｔｒｕｃｔｅｄ ｔｈｅ ｃｏｒｒｅｓｐｏｎｄｉｎｇ ｃｏｍｐｕｔｅｒ ｓｏｆｔｗａｒｅ Ｂｙ ｍｅａｎｓ ｏｆ ｎｕｍｅｒｉｃａｌ ｃｏｍｐｕｔａｔｉｏｎｓ，ｔｈｅ ｔｈｅｓｉｓ ｉｎｖｅｓｔｉｇａｔｅｄ ｔｈｅ ｓｔａｂｉｌｉｔｙ ａｎｄ ｔｈｅ ａｂｓｏｒｂｉｎｇ ｐｒｏｐｅｒｔｙ
时域有限差分法论文
时域有限差分法1 选题背景在多种可用的数值方法中，时域有限差分法(FDTD)是一种新近发展起来的可选方法。
1966年，K.S.Yee 首次提出电磁场数值计算的新方法—时域有限差分法(Finite Difference- Time Domain ，简称FDTD)。
经历了二十年的发展FDTD 法才逐渐走向成熟。
上世纪80年代后期以来FDTD 法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。
FDTD 法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将Maxwell 旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。
是电磁场和电磁波运动规律和运动过程的计算机模拟。
原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。
现在FDTD 法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域[1]。
经过了近四十年的发展，FDTD 法在计算方法和应用上取得了大量成果。
近几年来，讨论FDTD 法的深入发展和实际应用的文章几乎按指数增长，目前FDTD 法的许多重要问题得到了很好的解决，已经发展成为一种成熟的数值计算方法。
随着计算机数据处理性能的快速提高和计算机价格的下降，使得FDTD 法的应用范围越来越广，而FDTD 法本身在应用中又有新的发展.2 原理分析2.1 FDTD 的Yee 元胞E,H 场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理t t ∂∂=∂∂=⨯∇E D H ε t t ∂∂-=∂∂-=⨯∇HB E μ图1 Yee 模型如图1所示，Yee 单元有以下特点[2]：1）E 与H 分量在空间交叉放置，相互垂直；每一坐标平面上的E 分量四周由H 分量环绕，H 分量的四周由E 分量环绕；场分量均与坐标轴方向一致。
微带天线时域有限差分法分析
图 １天 线 输 入 阻 抗 随 馈 电 位 置 的变 化
（ ３ ） 易于集成 ，节省空间，降低成本
其存在的缺点主要有： （ １ ） 微带天线是 窄带器件，难于做到 宽频 带，难 于实
２ 微带印制天线的Ｆ Ｄ Ｔ Ｄ 建模
时域有 限差分 法尤其 适合计算 天线及 其他结 构 的宽
现 宽频带 跳频通信 。通 过增大介 质基板 的厚 度或者 是采 用 高介 电常数 的介质基板可适 当增宽天线 的带宽 。
较低 ，而且最大增益实 际上受限制 （ 约为２ ０ ｄ Ｂ ）。 （ ４ ） 介质基板 的厚度和材 质难于做到均匀分布 ，因此 相 同材料 制作 的微带天线 的性 能可能会有 较大 的性能差 距，因此微 带天线 的组 阵较为 困难 。 （ ５ ） 天线 的馈线 与辐 射体 之间存在距离 ，可能存在表
面波 。
变天线 的尺寸对 天线 性能 的影 响是 显而易见 的 ，所 以天
线 的小型化 技术式 当今天线 设计 的一个难 点 。微 带天线
１ ． ２ 微 带天线馈 电基本方式 总体而 言，微带天线有三种馈 电方式［ ｚ ］ ：（ １ ） 同轴馈
的 出现是小 型化技术 发展 的载体达 到共形等 特点 。通 过对微 带天线 的结 构进 行改造 ，可 以改变天线 的工作 电流分布 ，进而 改变
电流 的 关 系 如 下所 示 …：
天 线 的尺 寸 。本 文在 分 析 时 假 定天 线 长 为７ ０ ｍ ｍ ，宽 为
１ ０ ｍ ｍ ，印制在 介 电常数 为２ ． １ 厚为 １ Ｏ ｍ ｍ 的介 质板 上 。两 个矩 形槽 的作用 是在 原来天 线 的结 构上 ，额外 增加 了一 个Ｌ Ｃ 谐振 电路 结构 ，因此 ，天线工 作 出现 了两个 谐振 频 率 ， 当谐振 频 点接近 或无缝 连接 时 ，天线的带 宽就被 展 宽 了 。文 中 的天线采 用 同轴 馈 电方 式 ，馈 电点 位于 丽个 开槽 的中间 。金属 贴片 中心至两 个开槽的距离 为Ｐ 。改变 开槽 的长度 、宽度及ｐ 可 以改变天线的谐振频率，进而改 变 天线 的带宽。
时域有限差分法的图形处理单元的加速
时域有限差分法的图形处理单元的加速摘要：时域有限差分法，即FDTD（Finite Difference Time Domain），是计算电磁学的一种重要方法。
作为一种天然的并行算法，它的计算过程可以划分为多个同时进行相似计算的子计算。
这个方法主要是把麦克斯韦方程在时间上和空间上进行差分化，并且通过时间领域上的更新来模仿电磁场的变化来计算问题，因而有利于解决很多电磁场问题。
而图形处理单元即GPU（Graphic Processing Unit）相对于CPU的高性能计算速度以及NVIDA公司生产的GPU特有的高并行结构，为时域有限差分的加速提供了可能。
关键字：时域有限差分法；图形处理单元；麦克斯韦方程；并行算法1 FDTD的基本原理FDTD算法是1966年K.S.Yee发表在AP上的一篇论文建立起来的，后被称为Yee网格空间离散方式。
核心思想是把带时间变量的麦克斯韦旋度方程转化为差分形式，故中文称之为“时有限差分法”。
麦克斯韦方程如下：其中H是磁场强度，E是电场强度，D是电位移，B是磁通量密度。
上述两个矢量方程描述了麦克斯韦方程中磁场与电场复杂交错的关系。
由于在三维空间中每个矢量方程又可以分解为三个标量方程，因此该方程组可以化为6个标量方程如下所示：这样的话便将问题的几何空间离散为空间网格，电场和磁场的分量便被置于空间离散的网格点上。
而这是FDTD计算的前提。
我们再用差分近似替代麦克斯韦这6个标量方程中的时间和空间导数，构造一系列方程，均以前一时间步电磁场瞬时值来“预测”后一时间步电磁场的瞬时值，由此构造时间不断向前推近的算法，来模拟时域中的电磁场变化过程。
1966年，Yee首次给出了麦克斯韦旋度方程的一组差分形式，这组方程在空间和时间上一离散的形式给出，使用的是中心差分法。
对空间（X轴方向）的中心差分法离散公式如下：对Y轴，Z轴方向的中心差分离散以此类推。
对时间的中心差分离散公示如下：图1为Yee单元网格的结构：我们根据6个麦克斯韦标量方程，结合上述中心差分法，便可以得到三维问题下的FDTD更新方程如下：其余的五个方程也如法可以写出，因此任何时刻可一次算出一个点，并行算法可计算出多个点。
基于时域有限差分法对金属光栅异常透射研究
基于时域有限差分法对金属光栅异常透射研究基于时域有限差分法对金属光栅异常透射研究摘要：金属光栅异常透射是材料科学领域中具有重要实际应用价值的现象。
本文通过基于时域有限差分法的数值模拟，研究了金属光栅的异常透射效应。
通过建立合理的模型并选择合适的参数，得出了不同条件下的金属光栅异常透射的特性。
研究结果表明，金属光栅的异常透射效应与金属材料的参数、光波的入射角度以及光波频率密切相关。
这一研究对于充分发挥金属光栅的应用潜力有重要意义，有助于优化金属光栅在传感器、能源等领域的应用。
1. 引言金属光栅异常透射是指在光波穿过周期性金属结构时发生的非常规透射现象。
由于金属光栅具有微纳尺度的周期性结构，入射光波在金属光栅上发生多次反射和衍射，导致光波透射波矢的变化，从而产生不同于普通均匀介质的透射行为。
金属光栅异常透射现象在光学传感器、光电子学、纳米光子学等领域得到广泛应用。
2. 方法本文采用时域有限差分法（FDTD）作为数值模拟方法，对金属光栅的异常透射进行研究。
FDTD方法能够对光波在金属光栅中的传播情况进行准确描述，并可以得到透射光强的时域和频域特性。
通过建立合理的数值模型，考虑金属光栅的周期结构、金属材料的特性以及入射光波的参数等因素，进行数值模拟计算。
3. 结果与讨论我们在模拟中选择了不同的金属材料和光波参数，研究了金属光栅异常透射效应对这些因素的响应。
通过计算分析，我们得到了金属光栅异常透射的衍射效应和透射波矢变化规律等重要特性。
实验结果表明，金属光栅的异常透射效应与金属材料的电导率、透射光波的入射角度以及光波频率密切相关。
不同的金属材料和光波参数将导致不同的异常透射行为，这一研究结果对金属光栅的设计和应用具有重要指导意义。
4. 应用潜力与展望金属光栅异常透射现象在传感器、能源等领域有着重要的应用潜力。
通过深入研究金属光栅的异常透射特性，我们可以更好地利用这一现象，设计并优化金属光栅在传感器、能源收集等方面的应用。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
时域有限差分法1 选题背景在多种可用的数值方法中，时域有限差分法(FDTD)是一种新近发展起来的可选方法。
1966年，K.S.Yee 首次提出电磁场数值计算的新方法—时域有限差分法(Finite Difference- Time Domain ，简称FDTD)。
经历了二十年的发展FDTD 法才逐渐走向成熟。
上世纪80年代后期以来FDTD 法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。
FDTD 法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将Maxwell 旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。
是电磁场和电磁波运动规律和运动过程的计算机模拟。
原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。
现在FDTD 法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域[1]。
2 原理分析2.1 FDTD 的Yee 元胞E,H 场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理t t ∂∂=∂∂=⨯∇E D H ε t t ∂∂-=∂∂-=⨯∇HB E μ图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.2 Maxwell 方程FDTD 的差分格式麦克斯韦第一、二方程 ⎪⎪⎩⎪⎪⎨⎧⎪⎭⎫⎝⎛+∂∂-=⨯∇+∂∂=⨯∇m t t J B E J D H （1）式中，J 时电流密度，反映电损耗，mJ 是磁流密度，单位2m V /，反映磁损耗。
主要与上式对应。
各向同性介质中的本构关系：H JE J H B E D m mγγμε==== （2）其中m γ是磁阻率，计算磁损耗的。
以H E ,为变量，在直角坐标中，展开麦克斯韦第一、二方程，分别为⎪⎪⎪⎩⎪⎪⎪⎨⎧+∂∂=∂∂-∂∂+∂∂=∂∂-∂∂+∂∂=∂∂-∂∂zz x y y y zx x xy z Et E y H x H E t E x H z H E t E z H y H γεγεγε （3） ⎪⎪⎪⎩⎪⎪⎪⎨⎧-∂∂-=∂∂-∂∂-∂∂-=∂∂-∂∂-∂∂-=∂∂-∂∂zm z x y y m y zx x m xy z Ht H y E x E H t H x E z E H t H z E y E γμγμγμ （4） 令()t ,z ,y ,x f 代表H E,在直角坐标中的任何一个分量，离散符号取为()()()k ,j ,i ft n ,z k ,y j ,x i f t z y x f n=∆∆∆∆=,,, （5）()t ,z ,y ,x f 关于时间和空间的一阶偏导数取中心差分近似为()()[]()()[]()()[]()()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧⎥⎦⎤⎢⎣⎡-∆≈∂∂--+∆≈∂∂--+∆≈∂∂--+∆≈∂∂+∆=∆=∆=∆=k ,j ,i f k ,j ,i f t 1t f k ,j ,i f k ,j ,i f z 1z f k ,j ,i f k ,j ,i f y 1y f k ,j ,i f k ,j ,i f x 1x f 21-n 21n t n t 21n 21n z k z 21n21n yj y 21n21n xi x （6） 可以看出，每一节点上沿某一方向场分量的一阶偏微分可以用在该方向上相邻两点的一阶中心差商来描述，将式(1)用一阶中心差商方程取代，整理后便得到一阶差分方程，它具有二阶精度[3]。
Yee 元胞如图1所示，规定为1）剖分节点与场分量所在棱边中点不同，场分量的位置，即H E,节点是Yee 元胞节点的相对位置，不需要单独编码；2）当空间存在媒质分界面时，场量自动满足场的连续性条件，2t1t 2t 1t H H ,E E ==电磁分量的取样方式不仅符合法拉第电磁感应定律和安培环路定律的自然结构，也符合麦克斯韦方程的差分计算。
其次，时间步长可以取为电磁波传播一个空间步长所需时间的一半，因此E 与H 在时间顺序上交替抽样，时间间隔相差半个时间步长。
2.3 一维问题均匀平面波（TEM 波）是一维问题，设电磁波沿z 轴方向传播，则00==z z , H E ，场量和介质参数均与x ，y 无关，即0y ,0x =∂∂=∂∂，麦克斯韦方程为ym y xxxyH γt H μz E γE t E εz H +∂∂=∂∂-+∂∂=∂∂-（7）和xm x y yy x H γtH μz E γE tE εz H +∂∂=∂∂-+∂∂=∂∂- （8）旋转坐标轴后可以只保留一组公式[4]，设保留（7） Yee 元胞如图2所示E xH yz图2 一维Yee 元胞差分格式为()()()()()() k H k H z 1m CB k E m CA k E2121n y2121n y nx1n x⎪⎪⎭⎫ ⎝⎛--+∆-=+++ （9）()()()()()()[]k E 1k E z1m CQ k H m CP k H n x n x2121n y2121n y-+∆-+=+-+（10）如果介质无损耗，则0 ,0m ==γγ2.4 二维问题三维通常是散射问题，二维是TE 、TM 波问题，一维是TEM 波问题。
在二维场中，所有物理量与Z 坐标无关，既0z /=∂∂。
于是在TE 和TM 波的表达式分别为TE 波（0E z =） ⎪⎪⎪⎩⎪⎪⎪⎨⎧-∂∂-=∂∂-∂∂+∂∂=∂∂-+∂∂=∂∂H t H y E x E E t E x H E tE y H zm z x y y y zx x z γμγεγε （11）TM 波（0H z =） ⎪⎪⎪⎩⎪⎪⎪⎨⎧+∂∂=∂∂-∂∂-∂∂-=∂∂--∂∂-=∂∂zz x y y m y zxm x zEt E y H x H H t H x E H t H y E γεγμγμ （12）图3分别给出了TM 波和TE 波的Yee 元胞图图3 TM 波的Yee 元胞 图4 TE 波的Yee 元胞对于TE 波，只要令0=z E ，在z ∆上，yx H H , 不随z 变化，m 中去掉k 即可得到：()()()()()()j ,i H j ,i H y 1m CB j ,i E m CA j ,i E212121n z212121n z 21n x211n x⎪⎪⎭⎫ ⎝⎛-+-++∆++=++++ 式中：j ,i m 21+= （13）()()()()()()⎪⎪⎭⎫ ⎝⎛+--++∆-+=++++212121n z212121n z 21n y211n yj ,i H j ,i H x 1m CB j ,i E m CA j ,i E式中：21j ,i m += （14）()()()()()()()()⎥⎥⎦⎤⎢⎢⎣⎡∆+-++-∆+-++⋅-++=++-+yj ,i E 1j ,i E x j ,i E j ,1i E m CQj ,i H m CP j ,i H 21n x21n x21n y21n y212121n z212121n z式中，2121j ,i m ++= （15） 对TM 波，只要令0=z H ，在z ∆上，yx E E , 不随z 变化，m 中去掉k ，即可得到：()()()()()()[] j ,i E 1j ,i E y1m CQ j ,i H m CP j ,i H n z n z2121n x2121n x-+∆-+=+-+式中，21j ,i m += （16） ()()()()()()[]j ,i E j ,1i E x1m CQ j ,i H m CP j ,i H n z nz 2121n y2121n y-+∆++=+-+式中，j ,i m 21+= （17）()()()()()()()()⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡∆--+-∆--++=+++++yj ,i H j ,i H x j ,i H j ,i H m CBj ,i E m CA j ,i E 2121n x 2121n x 2121n y 2121n y n z 1n z式中：j ,i m = （18）为了编写统一的TE 和TM 波二维FDTD 程序，可将描述TE 波差分公式(13)～(15)中相应的标号整体移动1/2，即坐标（x,y ）分别沿x 和y 轴方向移动半个网格，并将离散时间也移动半个时间步长，式(13)～(15)可以重新写为()()()()()()[] j ,i H 1j ,i Hy 1m CB j ,i E m CA j ,i E n z n z2121n x2121n x-+∆++=+-+式中：21j ,i m += （19）()()()()()()[]j ,i H j ,1i H x1m CB j ,i E m CA j ,i E n z n z 2121n y2121n y-+∆-+=+-+式中：j ,i m 21+= （20）()()()()()()()()⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡∆--+-∆--+-=+++++yj ,i E j ,i E x j ,i E j ,i E m CQj ,i H m CP j ,i H 2121n x 2121n x 2121n y 2121n y n z 1n z式中，j ,i m += （21） 可以看出，TE 波的FDTD 公式(19)～(21)与TM 波的FDTD 公式(16)～(18)形式相同，给编程带来极大方便。