FDTD

FDTD（精选7篇）

FDTD 篇1

摘要：介绍了光子晶体波导的原理,然后运用时域有限差分方法(FDTD法)分析了二维光子晶体波导的传输特性,最后运用Matlab语言实现了二维光子晶体波导的仿真,使用的是TM模。

关键词：光子晶体,光子晶体波导,时域有限差分算法,光子带隙

0 引言

光子晶体的概念是由E.Yablonovitch[1]和S.John[2]于1987年几乎同时分别独立地提出的。自从光子晶体这个概念被提出以来,人们就在理论上和实验上对光子晶体进行了广泛的研究[3]。光子晶体是由不同介电系数的物质周期性排列而成的,在这种折射率呈周期性排布的电介材料中,某些波段的电磁波因周期性结构的强散射效应(strong scattering effect)将无法在电介材料中传播,因而形成光子带隙结构,即“光子禁带”。频率落在光子带隙内的电磁波不能在光子晶体中传播,因而能很好地控制光在介质中传播。

本文正是利用这一特性制成了光子晶体波导,它有效地为光通信系统中现存的问题提供了一个良好的解决途径,为将来超大规模全光或光电子集成回炉提供了一个很好的基础物质平台。

本文利用时域有限差分法(Finite Difference Time Domain ,简称FDTD 法)[4,6]分析了光子晶体波导的传输特性、模场分布还有透射率,时域有限差分法计算较为简便、通用性和适用性强、节约计算空间和存储空间。FDTD直接把含时间变量的Maxwell旋度方程在Yee氏网格空间中转换为差分方程,在时间和空间轴上逐步推进的求解,最终求出空间场的分布。

1 FDTD算法

FDTD算法: 把含时间变量的Maxwell旋度方程在Yee氏网格空间中转换为差分方程。在这种差分格式中每个网格点上的电场(或磁场)分量仅与它相邻的磁场(或电场)分量及上一时间步该点的场值有关。在每一时间步计算网格空间各点的电场和磁场分量,随着时间步的推进,即能直接模拟电磁波及其与物体的相互作用过程。

TM情形:直角坐标系下

undefined

undefined

将各式写成离散式:

undefined

2 介质的边界条件

Yee网格与垂直的X-Y平面如图1和图2所示。

最边缘的Hx,Hy,Ez由于没有更边缘的电磁场值可供计算,因此确定零为其默认值。

使用完美匹配层(Perfectly Matched Layer, PML)来作为吸收边界,完美匹配层设置的基本结构如图3所示, 在FDTD计算区域中,麦克斯韦方程以常规的FDTD方法求解。在计算域四周PML层,计算域中的散射体或者辐射源产生的外行波穿过与PML层的分界面,在PML层中被吸收。

整个网格分为九个区域,最中间是计算区域,区域1～8是APML层。其中区域5～8是APML的角区域,相应的电磁参数设置是(σx,σy,σmx,σmy)。区域1～2垂直x轴, 电磁参数设置是(σx,σy,0,0)。蛆域3～4垂直于y轴,电磁参数设置是(0,0,σmx,σmy)。

网格的划分过程中,时间步长取为undefined,其中c为自由空间光速。

3 仿真算法和结果分析

光子晶体参数:介质柱的晶格常数a=10.66665mm,半径R=2.13333mm(占空比f=0.2),其相对介电常数εa=11.7(GaAs)衬底相对介电常数为εb=1(空气)的正方晶格光子晶体。使用0.5mm×0.5mm的网格对21×21个晶包进行划分,四周用完美匹配层对光波进行吸收。首先对没有加入缺陷的光子晶体进行分析,激励源采用高斯脉冲,设置PML边界层数为12a,得到这种结构的光子晶体的带隙图, 得到光子带隙在0.27到0.42之间。

然后在光子晶体中引入线缺陷,构成光子晶体光波导, 线缺陷光子晶体光波导是在光子晶体的中间用空气取代一行介质柱。如图5所示。

在A点激励高斯脉冲,在B1点进行接收,通过对接收到的信号进行傅里叶变化后,比较接收端和输入端的光功率,对得到的图形进行归一化处理,得到如图6所示的透射系数。

从图中可以明显看到,原来光子带隙中的光可以较好地通过光子晶体光波导。

图7是其模场分布图。从图7中看到光子禁带的光频率只能在预先设计好的通道中通过,较好验证了光子晶体波导理论:当把缺陷引入光子晶体中,那么落在禁带频率的光波就被局限在这个缺陷中,如果引入的缺陷是线缺陷,那么就能将这些频率的光波从一个位置导向另一个位置,这就是一种全新的导波机制——光子带隙导波。

参考文献

[1]Yablonovitch E.Inhibited Spontaneous Emission in Solid-state Physicsand Electronics[J].Phys.Rev.Lett.,1987,58(20):2059-2061.

[2]John S.Strong Localization of Photons in Certain Disordered DielectricSuper Lattices[J].Phys.Rev.Lett.,1987,58(23):2486-2489.

[3]Jannopoulos J D,et al.Photonic crystals:modeling the filow of light[M].New York:Princeton university press,1995:12-13.

[4]葛德彪,闫玉波.电磁波时域有限差分方法[M].西安:西安电子科技大学出版社,2002.

[5]王长清,祝西里.FD-TD的基本原理[J].无线电电子学汇刊,1988(1-2):38-48.

[6]Yee K S.Numerical Solution of Initial Boundary Value Problems In-volving Maxwell’s Equations in Isotropic Media[J].IEEE.Trans.An-tennas Propagate.,1966,14(4):302-307.

[7]刘青,李承芳,张翔.二维光子晶体带隙结构的研究[J].激光杂志,2002.

FDTD 篇2

1 圆形槽波导单管座结构

1.1 管座位于腔体中心的输入阻抗

传统的波导振荡器, 一般管座位于波导腔体的中心。圆形槽波导中心单管座结构示意图如图1所示。偏置杆位于圆形槽波导的中央, 有源器件放在偏置杆的中间位置。

首先对该结构进行网格划分, 圆形槽波导的曲线导体边界的网格划分采用非均匀网格划分方法, 既可以描述局部细微结构, 又使运算量不至于很大。具体划分结构如图2所示。

图3为FDTD计算中使用的仿真模型, 圆形槽波导平板区开放边界上设置x方向的PML, 其参数为: (σx, σx*, 0, 0, 0, 0) , 用带左划线的阴影区表示;圆形槽波导z方向上设置z方向的PML, 其PML参数为: (0, 0, 0, 0, σz, σz*) , 用带右划线的阴影区表示;两种PML交界处的PML参数为: (σx, σx*, 0, 0, σz, σz*) , 用带花格线的阴影区表示;偏置杆位于圆形槽波导腔体的中心位置, 并且离开x方向的PML一定距离 (PML内边界不能与导体接触) 。

实际计算中采用的尺寸为:

中心槽区直径2a=14.9mm;

平板间距2c=9.7mm;

槽波导高度h=100mm;

x方向PML厚度h0=7.85mm;

槽波导长度l=50mm;

z方向PML长度l0=12.4mm;

管座直径3.5mm;

杆座距x方向PML内边界距离0.1mm。

考虑到数值色散问题, 圆形槽波导内部采用了较细致的网格划分, 中心槽区的网格大小大约是自由空间波长的1/20。

数值解的稳定性条件采用下式:

为了在边界上获得较理想的吸收效果, x方向的PML和z方向的PML均设置为18层。

首先设置强迫激励源。在中心管座处, 高斯脉冲源采用delta-gap模型, 即高斯脉冲源Ex在管座处均匀分布, 其表达式为

其中, 高斯脉冲源的有效频谱的中心频率f0=36GHz, 选择在有源器件的工作频率上。

中心管座处的电流i (nΔt) 可利用安培定律, 在每个时间步长对磁场进行环路积分得到。

对电压i (nΔt) 和电流i (nΔt) 分别进行离散傅立叶变换, 可得到管座结构的输入阻抗:

由此可得到管座输入阻抗曲线, 如图4所示。从曲线可以看出, 在激励源有效频谱中心f0=36GHz附近, 管座的输入电阻和输入电抗值最大。

1.2 管座偏离腔体中央的输入阻抗

考虑到圆形槽波导的大尺寸特点, 可用来制作功率合成器件, 在腔体中放置多个有源器件。这样, 管座位置就有相应的变化, 适当偏离腔体中心。下面就分析管座位置由中心向两边偏移时输入阻抗的变化情况。如图5所示。偏置杆位于圆形槽波导中央偏左位置, 有源器件放在偏置杆的中间。

参照前面的分析方法, 圆形槽波导横截面的划分同样采用非均匀网格划分方法, 其示意图见图6。FDTD计算中采用的仿真模型, 管座偏离中心距离为2.4mm, 其它设置与中心管座结构的仿真模型相同。

按照与中心管座结构输入阻抗的类似计算方法, 可得到偏离中心管座结构的输入阻抗。

根据上述方法, 选取不同的偏离位置进行计算, 再由插值方法得到了管座输入阻抗随偏离中心距离的变化曲线, 如图7所示。可以看出, 在中心处, 管座输入阻抗电阻部分最大, 随着偏离中心距离的增大, 管座输入阻抗电阻部分陡然变小, 电抗部分亦然。

2 结论

随着管座偏离中心位置的增大, 阻抗有相应的减小, 这意味着实现阻抗匹配将相对容易。另外考虑到功率合成器件中二极管对位置的排列, 启发我们研究单管管座偏离腔体中心的圆形槽波导振荡器的特性。考虑到器件的体积和腔体的大小以及工艺的难度, 适当选择管座位置偏离波导中心的距离, 可以在保证器件的设计和安装下, 最大程度减小外电路阻抗, 以利于阻抗匹配。

参考文献

[1]Hong-Sheng Yang, et al."Circular groove guide for short millimeter and sub-millimeter waves, "IEEE Trans.Microwave Theory Tech., 1995, 43:324-330.

[2]B.Denecker, F.Olyslager, L.Knockaert, and D.D.Zutter, "Generation of FDTD subsell equations by means of reduced order modling, "IEEE Trans., Antennas Propagat., 2003, 51:1806-1817.

FDTD 篇3

关键词：左手材料,Drude模型,时域有限差分法,电波传播

0 引 言

物理学中,介电常数ε和磁导率μ是描述介质中电磁场性质最基本的两个物理量。在已知的物质世界中,对于普通的电介质而言,介电常数ε和磁导率μ都为正值,电场、磁场与波矢三者构成右手螺旋关系,这样的物质被称为右手材料(Right-Handed Materials,RHM)。所谓的左手材料(Left-Handed Materials,LHM)是指介电常数ε和磁导率μ同时为负的介质材料,也常被称为双负介质(Double Negative Materials,DNG),其特点是电场、磁场与波矢三者构成左手螺旋关系。左手材料是近年来国际物理学和电磁学的一个研究热点,其概念最初由前苏联物理学家 Veselago 于1968年提出并做了大量的理论性研究,指出了左手材料具有诸如负折射效应、逆多普勒效应等许多奇异的电磁特性[1],但由于自然界中没能发现ε和μ同时为负数的介质材料存在,所以他的研究结果在很长一段时间一直没有得到实验验证,也没能激起人们更多的兴趣。1999年,英国皇家学院Pendry等人相续提出了用周期性排列的金属棒和开口金属谐振环可以在微波波段分别产生等效负介电常数和等效负磁导率的思路,并提出了左手材料具有“完美透镜”特性的概念[2]。2001年,美国加州大学圣迭哥分校物理学家Smith教授等人首次成功地通过人工方法构造出了这种自然界中并不存在的材料,并且利用此介质进行了电波传播实验,通过实验观察到了负折射等一系列左手材料中电波传播的特殊现象[3]。这些研究成果在国际上引起了很大的反响,激起了更多学者对左手材料在各个领域可能产生的应用前景进行了深入的思考和研究,而电磁波在该材料中的传播特性显然是研究的重要课题。

目前,物理光学方法、矩量法、时域有限差分法(FDTD)、高低频混合方法等各种数值分析方法纷纷被用来仿真和分析左手材料中电磁波的传播特性[4,5],其中,时域有限差分法(FDTD)特别是基于Drude模型的FDTD方法是比较方便和有效的一种[3,6]。但在以往的研究当中,出于典型和计算方便等原因,常常在 Drude模型中将左手材料的损耗因子设为0,这显然不完全符合实际情况。本文在建立有损耗左手材料Drude模型的基础上,推导了其FDTD算式,并在不同损耗因子情况下,对LHM模型中TM波传播的基本特性进行了数值仿真。结果表明,有损耗LHM内TM波传播的方向、相速度方向等特性与无损耗LHM情况一样,TM波在进入有损耗LHM后,电场幅度会衰减至边界值的$1/\sqrt{1-n}$。此外,有损耗LHM平板也能对点源信号进行聚焦,但其聚焦效果随着损耗因子、平板厚度及波源到平板的距离改变而变化。

1 有损耗LHM的Drude模型

由于左手材料具有负的折射率,必然存在色散与吸收,因此,在利用FDTD仿真分析左手材料时,为避免上述提到的数值发散现象,一种基于等离子体概念,能对左手材料的介电常数ε及磁导率μ进行间接设置的所谓Drude模型得到了广泛的采用。在Drude模型里,左手材料的相对介电常数和相对磁导率可表示为[4,7]:

$\begin{array}{l}\epsilon {}_{\text{r}}=\epsilon {}_{\infty }-\frac{\omega {}_{\text{e}}^2}{\omega {}^2+\Gamma {}_{\text{e}}\text{j}\omega }\\ \mu {}_{\text{r}}=\mu {}_{\infty }-\frac{\omega {}_m^2}{\omega {}^2+\Gamma {}_{\text{m}}\text{j}\omega }\end{array}(1)$

式中:ε∞为直流相对介电常数;ωe为电等离子体的频率;Γe为电碰撞系数;ωm磁等离子体的频率;Γm为磁碰撞系数;Γe和Γm代表的是材料的损耗因子,为得到理想无损耗的左手材料模型,只要令Γe=Γm=0即可。

利用式(1)可以方便求取左手材料的折射率。为了简化分析,令ε∞=μ∞=1,ωe=ωm=ωρ,Γe=Γm=Γ,此时折射率可表示为:

$n=\sqrt{\epsilon {}_{\text{r}}\mu {}_{\text{r}}}=1-\frac{\omega {}_{\rho }^2}{\omega {}^2+\Gamma {}^2}+\text{j}\frac{\Gamma \omega {}_{\rho }^2}{\omega {}^3+\omega \Gamma {}^2}$

式中:ω=2πf为入射波角频率。由于Γ代表左手材料的损耗因子,其值通常为10-8～10-9数量级[8],在本文所研究的30 GHz频率波段,显然有ω≫Г,由此进一步推导可知上式的虚部趋0,折射率只剩下实部,并可化简为:

$n(\omega )\approx 1-\omega {}_{\rho }^2/\omega {}^2(2)$

所以,尽管在Drude模型中,折射率n(ω)为一复数,但在一定的频率和损耗因子条件下,左手材料折射率n(ω)可近似为一实数,合理地对Drude模型中的各参数进行设置,可以得到折射率为不同负实数值的左手材料。

1.1 有耗LHM模型的FDTD算式

考虑本构关系:D=εrε0E,B=μrμ0H和极化电流密度J、极化磁流密度K,并采用变换∂/∂t↔-jω后,可得出左手材料广义Maxwell方程的时域形式:

基于简单和典型的考虑,本文主要分析和推导二维TM波情况。在此情况下,∂/∂z=0且Hz=Ex=Ey=0。构造Drude模型下二维TM波YEE元胞,用中心差分代替偏微分,将Durde模型εr和μr表达式代入并做整理,可得出二维TM波Drude色散媒质模型的FDTD算式。其中,Ez和Jz的表达式如下:

1.2 Drude模型的PML边界条件

在FDTD 算法中,为了在有限的时间和有限的计算机存储容量下模拟电磁波在无限大空间的传播,必须在FDTD模型中把网格截断,使网格空间成为有限。然而,为了减少误差在计算中,必须用某种合理的算法去模拟这种截断,即设置边界条件。目前,基于在截断边界处设置假想的可吸收入射波材料思路的FDTD边界条件,即完美匹配层(PML)边界条件,由于其计算的精确性得到了广泛的应用。在Drude模型下的PML层里,对于TM波,把电场Ez分解为两个子分量Ezx和Ezy,并且Ez=Ezx+Ezy,由PML层中TM波的Maxwell方程可得:

并且:

将以上频域方程转为时域方程,并进行差分离散,即得Drude模型的PML边界条件算式。

2 算例仿真结果与分析

2.1 损耗因子及折射率的影响

为简化分析,首先仿真一维情况。采用Matlab工具软件,用M语言编程,仿真一束频率f=30 GHz的平面电磁波垂直入射到无限大左手平板材料的情况。FDTD建模过程如下:一维空间被划分为x方向的200个网格,其中,在(100～150)Δx的区域放置无限大LHM平板。为保证迭代的收敛,时间步长Δt(单位:s)和空间步长Δx(单位:m)的选择必须满足Cournant稳定性条件,实际运算时取:

平面波激励源放置于LHM平板左侧20Δx处,采用PML边界条件,运算总时间步数取2 000。

为得到n=-1的典型左手材料,在Drude模型中选择ε∞=μ∞=1,等离子体角频率取$\omega {}_{\text{e}}=\omega {}_{\text{m}}=\omega {}_{\text{ρ}}=\sqrt{2}\omega$即可。图1、图2分别给出了在以上条件下,当LHM的损耗因子分别取Γe=Γm=10-6和Γe=Γm=10-9,程序运算到1 000步时Ez的瞬时值图。

从程序运行过程明显观察到:

(1) 在波源右侧的真空内,Ez的波峰随时间是向右推进的,表明真空中电磁波的相速度与能流密度传播方向一致。在有耗LHM平板内,Ez的波峰随时间向左推进,与传播方向相反,表明电磁波在有耗LHM内相速度为负值;

(2) 从图1和图2可以观察到,Ez波形在LHM和RHM两种材料分界面出现了尖点,这是因为Ez在边界上需要满足切向分量相等的边界条件,而电磁波在两种材料中的相速度方向相反造成的。这些结论与Smith等人在文献[3]中的分析及文献[4,9]等文献的结果相似。

对比图1和图2,随着损耗因子的增大,LHM内电场强度的振幅变小,且随传播距离的增加而衰减,这个结论与文献[5]的仿真结果正好相反。事实上,对于色散介质,其内部磁场能量密度为:

$w=w{}_{\text{e}}+w{}_{\text{m}}=\frac{\partial (\epsilon \omega )}{\partial \omega }E{}^2+\frac{\partial (\mu \omega )}{\partial \omega }\text{Η}\text{z}$

将Drude模型的关系式(1)、式(2)代入并考虑:

$\mu {}_{\infty }=\epsilon {}_{\infty }=1‚\Gamma {}_{\text{e}}=\Gamma {}_{\text{m}}=\Gamma ‚\omega {}_{\text{e}}=\omega {}_{\text{m}}=\omega {}_{\text{ρ}}=\sqrt{1-n}\omega$,可得:

而在RHM材料中,始终有:

显然,要保证能量守恒,LHM材料中的场值应该是RHM中场值的$1/\sqrt{1-n}$倍,因此,进入LHM的电波强度就会受到了抑制。

根据式(2),改变ωρ,还可得到不同折射率的LHM。图3和图4给出了固定损耗因子Γe=Γm=10-6而折射率分别为n=0和n=-6时,电场瞬时值的仿真图。

从图3看到,Ez在零折射率材料内几乎是一条直线,这表明电磁波在零折射率LHM中传播时,其相位是不发生变化的。图4则显示在n=-6的LHM内,其内部的场值是外部RHM场值的$1/\sqrt{7}$倍,这与式(3)推出的结论一致。

2.2 二维仿真

仿真一束频率f=30 GHz的波束垂直入射到有耗LHM平板的情况。二维空间被划分为600Δx×800Δy个YEE网格,取Δx=Δy= 10-4 m,Δt=1.67×10-13 s。激励源位于平板正上方[400,700]位置,LHM平板置于x=[0,800],y=[220,420]的网格区域内,其厚度与其上表面到波源的距离相等,均为200Δx=2λ0,边界上采用三层PML匹配层,厚度取10Δx,采用TM波多循环m-n-m脉冲[6]为激励源。对式(4)决定的FDTD算式用Matlab工具软件编程计算。取$\omega {}_{\text{e}}=\omega {}_{\text{m}}=\omega {}_{\text{ρ}}=\sqrt{2}\omega$,相应的折射率n=-1,加入耗散因子Γe=Γm=10-8 ,运算总时间步数为8 000。图5是运算时间步N=3 100时,电场Ez的伪彩色显示图。从该图明显看到:

(1) 电磁波的等相位面在LHM的内外是不同的,表明在LHM内电磁波相速度与传播方向相反,这与一维仿真结论一致;

(2) 在LHM内部,逐渐发散的电磁波又被重新聚集起来,最终在LHM内形成一个焦点,这个现象显然支持LHM的负折射理论;

(3) 在LHM另一侧的真空中,TM波再次出现了明显的聚焦现象。

以上仿真结果与J.B Pendry 等人在文献[2,3,6]中提到的左手材料具有“完美透镜”特性及能实现“二次聚焦”等结论相吻合。

进一步的仿真研究还表明,以上“二次聚焦”现象并不是在任何情况下都会产生的,损耗因子、LHM平板厚度以及折射率的改变,都会影响LHM对电磁波的汇聚效果。当保持激励源至LHM表面的距离不变,而LHM的厚度由2λ0减少至只有一个λ0时, LHM内的焦点将消失。而逐渐改变折射率,例如折射率由-1逐渐朝负方向变至-2时,在LHM内的焦点也由一个明显的焦点变成一个光斑,而光斑最后也将消失,但电磁波在LHM的轴线上仍有明显的聚焦效应。

3 结 语

在Drude模型下,采用时域有限差分法,对电磁波在穿过有损耗左手材料时的基本特性做了一定的仿真研究。一维数值仿真结果表明,电磁波在经过有损耗左手材料时,其相速度大小与光速相同,而传播方向则与电磁波传播的方向相反。同时,在边界处,LHM内的电场强度幅度是RHM内的$1/\sqrt{1-n}$倍,而TM波在零折射率材料中传播时,其相位是不发生变化的,随着损耗因子的增大,LHM材料内场的电场强度减弱。此外,二维数值仿真结果验证了平板左手材料具有类似普通右手材料凸透镜的聚焦效果,在特定条件下,平板左手材料能实现二次聚焦。这些仿真结果既支持了相关文献的理论和结果,同时也验证了方法的可行性,对左手材料的理解和进一步研究提供了一定的参考。

参考文献

[1]VESELAGO V G.The electrodynamics of substances with simultaneously negative values ofεandμ[J].Sov.Phys.Usp,1968,10(4):509-514.

[2]PENDRY J B.Negative refraction makes a perfect lens[J].Physical Review Letters,2000(85):3966-3969.

[3]SHELBY R A,SMITH D R,Nemat-Nasser.Microwave transmisson through a two-dimensional isotropic left-handed metamaterial[J].Appl.Phys.Lett.,2001(78):489-491.

[4]林振,梁昌洪.负媒质模型的时域有限差分法分析[J].强激光与粒子束,2006(4):615-617.

[5]林振,刘松华,代喜旺,等.有损耗左手材料的时域电磁分析[J].光电子激光,2007(7):813-815.

[6]ZIOLKOWSKI R W.Pulsed and CW gaussian beam interac-tions with double negative metamaterial slabs[J].Optics Express,2003,11(7):681-685.

[7]王芋,候德亭,冯宇,等.自由空间中左手材料传输特性的仿真研究[J].强激光与粒子束,2007(1):147-150.

[8]胡晶磊,汪蓉,周东山,等.左手材料与负折射[J].化学进展,2007(1):813-818.

FDTD 篇4

关键词：光子晶体波导,光子带隙,慢光,有限时域差分法

0 引言

随着集成电路在速度和效率上的发展, 其技术的提高受到了量子效应和电子本身相互作用的限制。光子技术则是突破这些限制的有效手段之一:光子响应时间快, 频谱宽;相互作用力低, 能耗小。因此, 如何控制光以实现光计算中信息处理能力的进一步增强是当今研究的热点之一。光子晶体因为其出色的对光的操控能力, 自被提出以来, 便获得了广泛的关注。

1987年, E·Yablonovitch[1]和S·John[1]在研究抑制自发辐射和光子局域态时分别独立提出了光子晶体的概念。所谓光子晶体, 是一种新型人造的介电常数呈周期性排布的介质结构。在光子晶体中, 由于介电常数空间上的周期性分布, 光波通过其中时, 光波的色散曲线将形成带状周期性结构, 带与带之间存在能隙, 也就是所谓的光子禁带。光子禁带是光子晶体最本质的特征, 与半导体中的电子能带结构类似, 这种相似性的存在使得固体物理中许多的概念都能够引用到光子晶体中, 比如布里渊区、倒格子、色散关系等等。

在文中, 利用FDTD方法分析一种新型的光子晶体慢光波导的色散关系。结构如下。首先分析了麦克斯韦方程的FDTD方法, 接下来利用得到的方程数值仿真新型空气型光子晶体慢光波导, 得到布里渊区边界处平坦的色散曲线, 对应为慢光模式。

1 FDTD

时域有限差分法 (FDTD) 由Yee在1966年提出[3], 直接求解麦克斯韦旋度方程。FDTD将空间和时间进行离散化处理, 将各离散点上的数值解来逼近连续场域内的真实解, 经过时间演算, 可以计算出电磁波随时间演化的规律, 给出精确的结构。是一种有效的求解电磁波方程的方法, 可适用于任何复杂的界面, 可以得到任何形状结构内部的电磁波。

在光子晶体中, 假设自由电荷ρ和电流$\overrightarrow{J}$都为零, 介质的参数不随时间变化并且各向同性, 那么麦克斯韦方程可表示为[4]:

$\begin{array}{l}\frac{\partial \overrightarrow{{\rm H}}}{\partial t}=-\frac{1}{\mu }[\nabla \times \overrightarrow{E}+\rho \overrightarrow{{\rm H}}]\\ \frac{\partial \overrightarrow{E}}{\partial t}=\frac{1}{\epsilon }[\nabla \times \overrightarrow{{\rm H}}-\sigma \overrightarrow{E}](1)\end{array}$

麦克斯韦旋度方程可以分解为下面六个耦合方程:

$\begin{array}{l}\frac{\partial {\rm H}{}_x}{\partial t}=\frac{1}{\mu }[\frac{\partial E{}_y}{\partial t}-\frac{\partial E{}_z}{\partial t}-\rho {\rm H}{}_x]\\ \frac{\partial {\rm H}{}_y}{\partial t}=\frac{1}{\mu }[\frac{\partial E{}_z}{\partial t}-\frac{\partial E{}_x}{\partial t}-\rho {\rm H}{}_y]\\ \frac{\partial {\rm H}{}_z}{\partial t}=\frac{1}{\mu }[\frac{\partial E{}_x}{\partial t}-\frac{\partial E{}_y}{\partial t}-\rho {\rm H}{}_z](2)\\ \\ \frac{\partial E{}_x}{\partial t}=\frac{1}{\epsilon }[\frac{\partial {\rm H}{}_z}{\partial t}-\frac{\partial {\rm H}{}_y}{\partial t}-\sigma E{}_x]\\ \frac{\partial E{}_y}{\partial t}=\frac{1}{\epsilon }[\frac{\partial {\rm H}{}_x}{\partial t}-\frac{\partial {\rm H}{}_z}{\partial t}-\sigma E{}_y]\\ \frac{\partial E{}_z}{\partial t}=\frac{1}{\epsilon }[\frac{\partial {\rm H}{}_y}{\partial t}-\frac{\partial {\rm H}{}_x}{\partial t}-\sigma E{}_z](3)\end{array}$

上面六个耦合偏微分方程是FDTD方法的基础, 分别对应TE模以及TM模。Yee在空间上建立了矩形差分网格 (如图1所示[5]) , 网格节点与一组对应的整数标号一一对应, 利用二阶精度的中心差分近似把旋度方程中的微分算符直接转换为差分形式。空间上使电磁分量交错放置, 在时间轴上电分量和磁分量相差半个时间步长。所采用的结构具有TE模式, 因此, 对式 (2) 进行了展开, 改写成差分方程:

$\begin{array}{l}{\rm H}{}_x^{n+\frac{1}{2}}(i,j+\frac{1}{2},k+\frac{1}{2})=\frac{1-\frac{\rho (i,j+\frac{1}{2},k+\frac{1}{2})△t}{2\mu (i,j+\frac{1}{2},k+\frac{1}{2})}}{1+\frac{\rho (i,j+\frac{1}{2},k+\frac{1}{2})△t}{2\mu (i,j+\frac{1}{2},k+\frac{1}{2})}}{\rm H}{}_x^{n-\frac{1}{2}}(i,j+\frac{1}{2},k+\frac{1}{2})+\\ \frac{△t}{\mu (i,j+\frac{1}{2},k+\frac{1}{2})}\cdot \frac{1}{1+\left\{\rho (i,j+\frac{1}{2},k+\frac{1}{2})△t/2\mu (i,j+\frac{1}{2},k+\frac{1}{2})\right\}}\times \\ \left\{\frac{E{}_y^n(i,j+\frac{1}{2},k+1)-E{}_y^n(i,j+\frac{1}{2},k)}{△z}+\frac{E{}_z^n(i‚j‚k+\frac{1}{2})-E{}_z^n(i,j+1‚k+\frac{1}{2})}{△y}\right\}(4)\end{array}$

在这里给出Hx分量的展开形式, 其余两个分量形式类似。从式 (4) 中可以看出, 每个网格点上的场分量的新值依赖于该点在前一时间步长时刻的值以及该点周围临近点上另一个场分量早半个步长时的值。因此任一时刻可以算出一个点, 利用并行算法可计算出周围的多个点。这样, 可以交替算出各个网格的电磁场, 最终得到想要时域数值结果。

在FDTD中, 时间间隔与空间间隔的选取需综合考虑:减小网格尺寸诚然可达到更高的精确度, 然而这也意味着计算存储量的大幅度增加。一般来说, 考虑到解的稳定性因素, 空间与时间离散间隔之间应当满足Courant稳定性条件:

$c△t\le 1/\sqrt{\frac{1}{(△x){}^2}+\frac{1}{(△y){}^2}+\frac{1}{(△z){}^2}}(5)$

其中, c为真空中的光速。

另外, 在数值模拟电磁场的辐射和散射问题中, 边界是开放的, 电磁场占据了无限大的空间;而计算机的内存是有限的, 所以只能模拟有限的空间:时域有限差分的网格在某处将被截断。这样就要求在边界处不能引起波的明显反射, 才能对向外传播的波而言就如在无限大的空间上传播。根据所计算的结构特点, 可以采用布洛赫周期边界以及吸收边界等方法。在所设计的光子晶体波导结构中一般采用后者。所谓吸收边界, 即波传到边界处被吸收而不产生发射。目前提出的吸收边界条件有很多, 在计算中, 采用了1994年由Berenger提出的完全匹配层 (PML, perfectly matched layers) [6], 具体方法及原理可参考文献[4]中所介绍。完全匹配层是指将电磁场分量在吸收边界区分裂, 并分别对各个分裂的场分量赋以不同的损耗, 这样在外边界处得到了一种非物理的吸收材料。

2 慢光光子晶体的FDTD仿真

所采用的结构如图2中所示, 此结构为三角晶格的二维光子晶体, 首先在中心去除一排空气孔形成普通的W1型波导, 再在波导正中心插入空气槽形成空气型波导。假设光子晶体的晶格常数为Λ, 取空气孔半径r/Λ=0.3, 中心空气槽宽度为0.2Λ。经过PWE方法已经算得, 此结构的带隙范围为0.202-0.275 (2πc/α) 。

取20个单位长度的光子晶体慢光波导进行分割计算, 采用晶胞单元为x方向上为20一个晶格常数, y方向上为9个晶格常数, 将晶胞划为分640×288个格子进行计算, 在结构四周采用厚度为Λ的完美匹配层, 在波导左边处 (图中红色线段表示) 设置高斯脉冲源。鉴于所采用的结构中存在TE模式, 因此设置高斯脉冲为TE模, 中心频率假设为0.25c/Λ。得到的色散关系曲线为图3所示。图中椭圆形标示的为空气芯所对应的模式, 对比普通W1型波导的色散曲线可以看出, 此能带为典型光在低介电常数介质中传播时所对应的色散曲线, 其斜率为正值。并且此时模式落在带隙范围内, 其接近于布里渊区边界的较宽范围内相对比较平坦。这意味着在此处可以得到慢光。

3 结束语

光子晶体波导是一种有效的获得慢光的结构。在本文中, 设计了一种新型的光子晶体结构, 引入FDTD方法对其色散关系进行计算。得到的结果与之前所采用频域的平面波展开法得到的结果非常一致。在布里渊区边界处得到了较宽范围内的平坦色散曲线, 对应慢光模式。

参考文献

[1]Yablonovitch E.Inhibited spontaneous emission in solid-state phys-ics and electronics[J].Phys.Rev.Lett, 1987, 58:2059-2062.

[2]John S.Strong localization of photons in certain disordered dielectricsuper lattices[J].Phys.Rev.Lett, 1987, 58:2486-2489.

[3]Kane Yee.Numerical solution of initial boundary value problems in-volving Maxwell's equations in isotropic media[J].Antennas andPropagation, IEEE Transactions, 1966, 14:302-307.

[4]葛德彪.电磁波时域有限差分方法[M].西安:西安电子科技大学出版社, 2002.

[5]Steven G.Johnson, Yee lattice.[EB/OL]http://en.wikipedia.org/wiki/File:Yee-cube.svg, 2005.

FDTD 篇5

文中研究THz波段薄层金属平板的屏蔽特性。在THz波段,室温下的金属电导率将采用经典弛豫效应模型[4]。在FDTD方法仿真中,采用Z变换实现金属色散媒质本构关系从频域到时域的转化,有效实现了FDTD方法中场方程的迭代更新。最后通过数值仿真,分析了不同金属材料对THz波的屏蔽特性,同时基于屏蔽效应大小与金属平板厚度的关系,此方法有望用于小于微米量级金属平层厚度的测量。

1 THz波段金属电导率模型

在微波频率,金属导体往往用复数等效介电常数表征[4]

$\epsilon {}_{eq}=\epsilon {}_0\epsilon {}_r=\epsilon {}_0\left(1+\frac{\sigma }{\text{j}\omega \epsilon {}_0}\right)$ (1)

其中,ε0为自由空间介电常数;σ为金属电导率,在微波频率下为常数。在微波频率传导电流远大于位移电流(σ≫ωε0),因此可以忽略εr的实部。然而当信号工作频率达到THz波段时,金属电导率色散效应趋于明显,需要用更精确的模型描述金属电导率。

在THz波段,利用金属导体的自由电子气近似模型,可得到金属电导率的经典弛豫效应模型[5]

$\sigma =\frac{\sigma {}_0}{1+\text{j}\omega \tau }$ (2)

其中,σ0为金属固有直流电导率;τ为平均电子碰撞时间。将式(2)代入式(1),得到金属等效介电常数模型

$\epsilon {}_{eq}=\epsilon {}_0\epsilon {}_r=\epsilon {}_0\left[1+\frac{\epsilon {}_0/\epsilon {}_0}{\text{j}\omega +(\text{j}\omega ){}^2\tau }\right]$ (3)

式(3)给出的金属等效介电常数即熟知的Drude模型,当ωτ≪1时,式(3)将退化为微波波段下的金属等效介电常数模型式(1)。在THz波段式(3)的Drude模型可较好地表征金属导体的等效介电参数。室温(20 ℃)下金属Drude模型中的参数可由实验测量得到,表1给出了银、铜、金、铝4种金属的固有直流电导率σ0和平均电子碰撞时间τ的值[6]。

2 FDTD方法处理一维金属平板模型

一维金属平板模型如图1所示,在xy平面为无限大,在z方向厚度为T,位于0≤z≤T区域,这样图1所示的区域1、2、3分别为空气、金属平板和空气3种媒质。假设THz脉冲沿z方向自空气垂直入射到金属平板中,电场分量和磁场分量分别在x方向和y方向。金属平板中,一维情况的Maxwell方程如下

$\frac{\partial E{}_x}{\partial z}=-\mu {}_0\frac{\partial {\rm H}{}_y}{\partial t}$ (4)

$\frac{\partial {\rm H}{}_y}{\partial z}=-\frac{\partial D{}_x}{\partial t}$ (5)

其中,μ0为真空磁导率;Dx为电位移矢量分量。对金属导体考虑其等效介电常数的色散特性时,频域电位移矢量和电场强度的本构关系为

Dx(ω)=ε0εr(ω)Ex(ω) (6)

利用FDTD方法对式(4)和式(5)在一维空间差分离散[7]

$D{}_x^{n+1}(k)=D{}_x^n(k)-\frac{\text{Δ}t}{\text{Δ}z}[{\rm H}{}_y^{n+1/2}(k+1/2)-{\rm H}{}_y^{n+1/2}(k-1/2)]$ (7)

${\rm H}{}_y^{n+1/2}(k+1/2)={\rm H}{}_y^{n-1/2}(k+1/2)-\frac{\text{Δ}t}{\mu {}_0\text{Δ}z}[E{}_x^n(k+1)-E{}_x^n(k)](8)$

对频域本构关系式(6),相对介电常数εr(ω)可写成以jω为自变量的有理分式形式[8,9]

$\epsilon {}_r(\omega )=\frac{{\displaystyle \sum _{n=0}^{{\rm N}}p}{}_n(\text{j}\omega ){}^n}{{\displaystyle \sum _{n=0}^{{\rm N}}q}{}_n(\text{j}\omega ){}^n}$ (9)

由微分的时频域相互关系jω←→∂/∂t,可将式(6)的频域本构关系转化到时域

$D{}_x(t)=\epsilon {}_0\epsilon {}_r\left(\frac{\partial }{\partial t}\right)E{}_x(t)$ (10)

将式(9)代入式(10)

$\left[{\displaystyle \sum _{n=0}^{{\rm N}}q}{}_n\left(\frac{\partial }{\partial t}\right){}^n\right]D{}_x=\epsilon {}_0\left[{\displaystyle \sum _{n=0}^{{\rm N}}p}{}_n\left(\frac{\partial }{\partial t}\right){}^n\right]E{}_x$ (11)

在时域要实现式(11)场分量时间微分的差分离散。引入函数$y(t)=\frac{\partial f(t)}{\partial t}$在(n+1/2)Δt时刻进行中心差分近似

$\frac{(y{}^{n+1}+y{}^n)}{2}\approx \frac{(f{}^{n+1}-f{}^n)}{\text{Δ}t}$ (12)

引进移位算子Zt满足Ztgn=gn+1,g=y,f,代入式(12),有

$y{}^n\approx \left(\frac{2}{\text{Δ}t}\frac{{\rm Z}{}_t-1}{{\rm Z}{}_t+1}\right)f{}^n$ (13)

对比$y(t)=\frac{\partial f(t)}{\partial t}$和式(13),时域微分算子可用下面离散双线性变换近似

$\frac{\partial }{\partial t}\to \frac{2}{\text{Δ}t}\frac{{\rm Z}{}_t-1}{{\rm Z}{}_t+1}$ (14)

将式(14)代入式(11),可得到电位移矢量分量Dx和电场分量Ex在时域满足的方程

$\begin{array}{l}\left[{\displaystyle \sum _{l=0}^{{\rm N}}q}{}_l\left(\frac{2}{\text{Δ}t}\right){}^l({\rm Z}{}_t+1){}^{{\rm N}-1}({\rm Z}{}_t-1){}^l\right]D{}_x^n(k)=\epsilon {}_0\\ \left[{\displaystyle \sum _{l=0}^{{\rm N}}p}{}_l\left(\frac{2}{\text{Δ}t}\right){}^l({\rm Z}{}_t+1){}^{{\rm N}-1}({\rm Z}{}_t-1){}^l\right]E{}_x^n(k)(15)\end{array}$

将金属的等效介电常数式(3)与式(9)对比,可知N=2,基于金属的等效介电参数模型式(3)及式(15),可得金属材料中电场E${}_x^{n+1}$(k)的更新方程

$\begin{array}{l}E{}_x^{n+1}(k)=\frac{1}{b{}_0}[\frac{a{}_0}{\epsilon {}_0}D{}_x^{n+1}(k)+\frac{a{}_1}{\epsilon {}_0}D{}_x^n(k)+\frac{a{}_2}{\epsilon {}_0}D{}_x^{n-1}(k)-\\ \begin{array}{l}\\ \end{array}b{}_{1}E{}_x^n(k)-b{}_{2}E{}_x^{n-1}(k)](16)\end{array}$

其中,$a{}_0=q{}_0+q{}_1\frac{2}{\text{Δ}t}+q{}_2\left(\frac{2}{\text{Δ}t}\right){}^2,a{}_1=2q{}_0-q{}_2\left(\frac{2}{\text{Δ}t}\right){}^2,a{}_2=q{}_0-q{}_1\frac{2}{\text{Δ}t}+q{}_2\left(\frac{2}{\text{Δ}t}\right){}^2$;$b{}_0=p{}_0+p{}_1\frac{2}{\text{Δ}t}+p{}_2\left(\frac{2}{\text{Δ}t}\right){}^2,b{}_1=2p{}_0-p{}_2\left(\frac{2}{\text{Δ}t}\right){}^2,b{}_2=p{}_0-p{}_1\frac{2}{\text{Δ}t}+p{}_2\left(\frac{2}{\text{Δ}t}\right){}^2$;$p{}_0=\frac{\sigma {}_0}{\epsilon {}_0},p{}_1=1,p{}_2=\tau$;q0=0,q1=1,q2=τ。

3 THz波段金属平板的屏蔽特性

分析THz波段金属平板的透射屏蔽特性,首先需要计算金属平板对THz波的透射系数

${\rm T}=\left|\frac{\overline{E}{}_t}{\overline{E}{}_i}\right|$ (17)

其中,$\overline{E}{}_i$为金属屏蔽物不存在时的入射场强;$\overline{E}{}_t$为存在金属屏蔽物时的透射场强。计算取距离金属平板右侧3个剖分网格处为透射系数观察点。另外计算透射系数时需将时域入射场强和透射场强利用傅里叶变换变换到频域。

屏蔽效应(Shielding Effects,SE)表征了金属导体对电磁波的传输隔离特性,其定义为[1]

$\text{S}\text{E}=20\log {}_{10}\left|\frac{\overline{E}{}_i}{\overline{E}{}_t}\right|$ (18)

由式(18)可见由透射系数T易于得到金属平板的屏蔽效应

SE=-20lgT (19)

3.1 不同金属材料的屏蔽特性

首先分析不同金属材料对THz波的透射屏蔽特性。入射波为微分高斯脉冲

$E{}_{xi}(t)=\left(\frac{t-t{}_0}{\tau }\right)\exp \left[-\frac{4\text{π}(t-t{}_0){}^2}{\tau {}^2}\right]$ (20)

脉冲宽度τ=170 fs,t0=260 fs。FDTD空间剖分网格Δz=37.5 nm,时间步长Δt=Δz/2c=0.062 5 fs。金属平板厚度取T=5Δz=187.5 nm,采用一阶Mur吸收边界条件,假设左右两侧边界剖分网格角标为1和K时,其边界更新方程分别为E${}_x^{n+2}$ (1)=E${}_x^n$ (2),E${}_x^{n+2}$(K)=E${}_x^n$(K-1)。

利用表1给出的金属材料参数,图2给出了不同材料金属平板的屏蔽效应。由图2可见,不同金属的屏蔽效应均随频率增大而增大,且具有周期振荡性。由于不同金属具有不同的固有直流电导率和平均电子碰撞时间,在频率较低时(<2 THz),其屏蔽主要受固有直流电导率的影响,银的固有直流电导率最大,其屏蔽效应最大,相应铝的屏蔽效应最小;随着频率升高,金属的屏蔽效应受平均电子碰撞时间的影响越来越显著,铝的平均电子碰撞时间最短,THz波在铝中的衰减随频率增大趋势较快,因此铝的屏蔽能力随频率增大较其他3种金属显著。4种金属在187.5 nm的厚度下屏蔽效应都>60 dB,可见金属平板对THz波的屏蔽能力很强,其主要能量被金属反射。

3.2 金属平板厚度对屏蔽效应的影响

取金属材料为金,分析金属平板不同厚度时对THz波的屏蔽效应。金属平板厚度取75～337.5 nm。图3给出了金属平板厚度分别取T=75 nm,150 nm,225 nm,300 nm时金属平板的屏蔽效应随频率的变化,由图3可见随着金属平板厚度增大,同一频率下的屏蔽效应逐渐增大。基于不同厚度金属平板对THz波屏蔽效应的大小,图4分别取固定频率f=5.05 THz(对应局部屏蔽极大值)和f=6.05 THz(对应局部屏蔽极小值),给出了金属屏蔽效应随金属平板厚度的变化曲线。由图可见,随着金属平板厚度的增大,金属平板对THz的屏蔽衰减分贝数线性增大。基于此特征,可考虑通过测定固定频率下金属平板的屏蔽效应大小来测量金属平板的厚度,这在测量小于微米量级的薄膜金属厚度中具有实际工程意义。

4 结束语

在分析一维金属平板的透射屏蔽效应中,考虑THz波段金属电导率的频率色散特征,在FDTD方法中利用Z变换实现频域同构关系的时域转化,使编程易于实现。仿真结果表明,对不同金属材料,其屏蔽效应均随频率增大且周期振荡,屏蔽大小受金属的固有直流电导率和平均电子碰撞时间的影响。同一频率下屏蔽效应随金属平板厚度增大线性增大,利用此特征可用于测量厚度小于μm级金属平板的厚度。

摘要：研究了室温下金属平板对THz波的透射屏蔽特性。考虑到THz波段金属电导率的频率色散特征,基于Z变换方法处理频域色散媒质本构关系到时域的转化,使FDTD仿真易于编程实现。仿真表明,对不同金属材料,其屏蔽效应随频率呈现出周期振荡,且金属平板厚度在百nm量级时对THz波的屏蔽效应>60 dB。对不同厚度的金属平板,其屏蔽效应以dB为单位随厚度线性增大,有望利用此特性实现厚度μ<m级薄层金属的厚度测量。

关键词：THz波,金属平板,色散模型,FDTD方法,屏蔽

参考文献

[1]LUCYSZYN S,ZHOU Y.Charactering room temperature THzmetal shielding using the engineering approach[J].Progressin Electromagnteics Research,2010,103(6):17-31.

[2]ANA V A,DAWOOD M,FALCONE F,et al.Analysis of timeand spatial evolution of Brillouin precursors through metallicmedia in THz band[C].Proceedings of the 5th EuropeanConference on Antennas and Propagation(EUCAP),2011:1608-1611.

[3]MAIK S,CHRISTIAN J,MARTIN K.Analyzing sub-100-umsamples with transmission terahertz time domain spectroscopy[J].Optics Communications,2009,282(6):1304-1306.

[4]CARLOS A L,JORGE A R,JOSE R M,et al.Field propaga-tion in circular hollow waveguides with non-ideal metallicconductors from microwaves to terahertz frequencies[J].IEEE Transactions on Microwave Theory And Techniques,2011,59(12):3013-3022.

[5]ZHOU Y,LUCYSZYN S.HFSSTM modeling anomalies withTHz metal-pipe rectangular waveguide structures at roomtemperature[J].Piers Online,2009,5(3):201-211.

[6]LUCYSZYN S.Investigation of anomalous room temperatureconduction losses in normal metals at terahertz frequencies[J].IEE Proc.Microwaves,Antennas and Propagation,2004,151(4):321-329.

[7]葛德彪,闫玉波.电磁场时域有限差分方法[M].西安:西安电子科技大学出版社,2005.

[8]葛德彪,吴跃丽,朱湘琴.等离子体散射FDTD分析的移位算子方法[J].电波科学学报,2003,18(4):359-362.

FDTD 篇6

目前, 并行算法的实现主要基于两种标准: (1) 基于分布式内存模型的MPI标准[3]; (2) 基于共享式内存模型的Open MP标准[4]。MPI标准下, 各个处理单元具有各自独立的局部储存器, 由于不存在公共可用的存储单元, 各处理器之间通过消息传递来交换信息, 控制并行任务的执行。Open MP标准下, 各个处理单元通过对共享内存的访问来实现信息的交换、任务的并行。

本文先以三维非磁化等离子体为例, 给出了处理色散介质的线性递归卷积 (PLRC) 方法。之后介绍了Open MP方法以及在色散介质时域有限差分方法求解中的具体实现步骤。最后计算了经过Open MP方法并行加速后, 在5~25 GHz频率范围内的电磁波对非磁化等离子体的透射衰减, 并对比了使用Open MP方法加速前后的计算时间。通过实验验证了该并行算法的有效性。

1 色散介质FDTD方法

1.1 色散介质基本模型

各向同性介质的麦克斯韦旋度方程表达式为[5]

对于色散介质, 其本构参数随时间变化而变化, 这里假定介质的磁本构参数与频率无关, 只有介电系数与频率有关。给出频域的本构关系式

本文研究的非磁化等离子体属于Drude模型, 其形式为

1.2 PLRC-FDTD方法

本文采用分段线性递归卷积时域有限差分方法 (PLRC-FDTD) 。与传统的线性递归卷积时域有限差分 (RC-FDTD) 相比, PLRC-FDTD假设每个网格的场值是线性变化的, 这使得PLRC-FDTD的计算精度相对于RC-FDTD来说有了较大改善。在PLRC-FDTD中, 卷积项中的电场场值在每一时间步长上是线性变化的。这里给出Ei+1、Ei、t的函数

将时间离散, 结合式 (4) 代入式 (3) , 得

以上是非磁化等离子体中电场的迭代方程, 磁场的迭代方程与传统FDTD方法一致。

2 Open MP并行设计

2.1 Open MP概述

Open MP本身并不是一种独立的语言, 而是对常用科学计算语言 (C/C++、Fortran) 的一种补充和扩展。Open MP通过编译指导指令和应用程序编程接口 (API) 来实现共享内存的多线程并行计算。

2.2 Open MP并行模型

Open MP的执行模型采用fork-join的形式, 其中fork创建新线程或者唤醒已有的线程;join即多线程的会合。fork-join执行模型在刚开始时, 只有一个被称为“主线程”的线程存在。在程序运行中, 当需要进行并行计算时, 则会调用新的线程去执行并行任务。在执行并行计算时, 主线程与调用的新线程同时工作。当并行任务结束时, 被主线程调用的新线程退出, 主线程依旧单独工作。因而, 在成对的fork和join之间的区域就是并行域, 放置需要并行执行的代码。

Open MP中常用的两种并行策略[6]:数据并行和功能并行。数据并行常使用for/parallel for指令制导, 功能并行常使用sections/parallel sections指令制导。

FDTD的迭代方式是蛙跳式 (forg-leap) En→Hn→En+1→Hn+1, 在每一时间步内, 电场的迭代与磁场的迭代相对独立, 即在同一时间步内, 只更新电场、磁场中的一项, 另一项保持不变。这样的迭代方式适合与并行计算。

在此选用数据并行的策略, 采用for与parallel for指令混合制导。

3 并行FDTD性能分析

3.1 数值算例

首先, 验证Open MP并行算法应用到色散介质FDTD方法时的准确性。这里以电磁波在三维分层非磁化等离子体中传播为例。入射波为平面波, 频率范围在5~25 GHz, 空间步长Δ=0.5 mm, 时间步长Δt=0.83 ps, 由连接边界条件引入。设置分层非磁化等离子体, 尺寸为5 cm×5 cm×5 cm, 沿入射方向分为3层。第1层厚度为1.25 cm, 等离子体参数为ωp=17.16 Grad/s, υ=6 GHz;第2层厚度为2.50 cm, 等离子体参数为ωp=25.74 Grad/s, υ=6 GHz;第3层厚度为1.25 cm, 等离子体参数为ωp=17.16 Grad/s, υ=6 GHz。在等离子体背对来波方向面的中心, 设置一个观测点, 对记录的透射波信号作傅里叶变换, 求得5~25 GHz频带内相对于入射信号振幅的衰减情况, 并将Open MP并行计算的结果同CST的仿真结果进行对比, 结果如图2所示。由对比可知, 加入Open MP并行算法后的计算结果与实际情况相符, Open MP并行求解色散介质的方法具有较高的数值精度。

3.2 并行性能分析

首先, 给出两个衡量并行性能的指标, 加速比和效率[7]。

加速比:单个线程运行程序花费的时间与n个线程参与运行程序花费时间的比值;效率:加速比与n的比值。

为了计算Open MP并行算法对于色散介质FDTD方法的并行性能, 这里给出了3组计算区域大小不同的模型, 计算结果收敛所需的迭代运行的次数分别为500次、1 000次、2 000次, 分别计算3组模型的加速比和效率, 3组模型在并行计算是均占用5个线程, 计算结果如表1所示。

4 结束语

本文利用Open MP编程, 实现了PLRC-FDTD方法的并行化计算, 并行计算的结果与CST仿真结果吻合, 具有较高的数值精度, 以及加速比和效率。因此, 大幅提高了FDTD方法针对色散介质的求解速度。

摘要：分析了色散介质时域有限差分的模型, 并针对非磁化等离子体给出了的分层线性递归卷积算法。介绍了Open MP并行设计的基本模型, 并将其应用于非磁化等离子体的计算当中。最终通过验证非磁化等离子体的透射电磁波, 验证了将Open MP并行设计应用于色散介质中的准确性, 同时依据计算区域大小不同的3个算例验证了该算法具有较高的并行性能。

关键词：时域有限差分,色散介质,OpenMP,并行程序

参考文献

[1]Kelly.Piecewise linear recursive convolution for dispersive media using FDTD[J].IEEE Transactions on Antennas and Propagation, 1996 (6) :792-796.

[2]Taflove.Computational electrodynamics:the finite-difference time-domain method[M].Boston:Artech House, 2000.

[3]莫则尧, 袁国兴.消息传递并行编程环境MPI[M].北京:科学出版社, 2011.

[4]Anon.Open MP application program interface, version2.5[EB/OL]. (2005-05-20) [2015-02-10]http://www.openmp.org.

[5]葛德彪.电磁波时域有限差分方法[M].西安:西安电子科技大学出版社, 2011.

[6]李正浩, 周俊, 遛达刚.基于Open MP的电磁场FDTD并行程序性能分析[J].现代电子技术, 2008, 31 (5) :135-136.

FDTD 篇7

一、一维FDTD差分方程表达形式

一维波动方程, 在一维无限长介质中, 其解为u=f (x-vt) , 是以速度v向x轴正向传播的行波。用差分形式, 把波动方程进行改写:

二、初始、边界条件的确定以及编程计算的结果

边界条件:u (t, 1) =0, u (t, Nx) =0, 对t=1到Nt。

初始条件:由于在差分形式的递推公式 (1) 中, 要计算t=3的值, 需要知道t=1, 2时的值。本文规定初始条件为u (t, i) =0, 对t=1, 2.

下面是两组不同的参数, 用Matlab编程计算的结果:

第1组:Lx=1×10-4, Δx=1×10-6, Nx=100, Nt=100, v=c=3×108,

信号源, 观察点在i=80, 也就是x=80Δx处。

第2组:Lx=1×10-4, Δx=1×10-6, Nx=100, Nt=1000, v=c=3×108,

信号源, 观察点在i=80, 也就是x=80Δx处。参数中时间长度Nt与第1组不同, 其他参数相同。

第3组:Lx=1×10-4, Δx=1×10-7, Nx=1000, Nt=3000, v=c=3×108,

信号源, 观察点在i=500, 也就是x=500Δx处。

从图中可以看出:

t=1到600, 在观察点无信号, 信号还没有到达;

t=600到900, 是过渡阶段, 与数值化的计算有关, 在观察点信号尚不稳定;

t=900到2400, 观察点信号与用连续方程得到的结果一致;

t=2400以后, 观察点的信号的幅度比源的信号还大, 是因为文中取一维两端边界为零, 反射波到达观察点, 与正向波形成干涉。

用快速傅里叶变换得出的频谱透射率, 由于已经包含了开始的过渡阶段和以后的反射波的干涉的影响, 在0.8左右, 与波的无衰减传播结果接近。

三、结论

在一维的FDTD计算中, 直接使用整数时刻和整数空间坐标, 而不需要使用半整数时间和半整数空间坐标, 简化了坐标的标记形式;时间步长取, 时间总长Nt的取值, 要求在这个时间信号源的信号可以到达观察点, 但又没有界面反射波到达观察点, 这样的参数比较理想, 可以得到与连续方程相同的结果。

摘要：在用一维时域有限差分方法计算波的传播时, 为了防止出现发散的结果, 要求时间步长Δt<Δx/v, 但具体取多大的值, 没有定论。通过编程试验, 时间步长取Δt<0.75Δx/c, 而时间步数Nt大于波从波源传到观察点所需要的时间, 又小于反射波到达观察点的时间。而且计算的时间格点和位置格点全部都可以取整数。

关键词：FDTD,有限时域差分方法,参数

参考文献

[1]葛德彪, 阎玉波.电磁波时域有限差分方法[M].第2版.西安:西安电子科技大学出版社, 2005.

[2]曾辉, 杨亚培.FDTD法与平面波展开法在光子禁带计算中的差异分析[J].电子科技大学学报, 2005, 34 (6) :901-904.

[3]张国华, 袁乃昌, 付云起.FDTD方法分析光子带隙能带结构[J].微波学报, 2001, 17 (4) :14-17.

[4]张亮, 寇晓艳.FDTD的Matlab语言的实现[J].延安大学学报:自然科学版, 2009, 28 (2) :57-59.