杨利霞 许红蕾 孙 栋 王洪金

(江苏大学计算机科学与通信工程学院通信工程系,江苏 镇江 212013)

引 言

时域有限差分(Finite-Difference Time-Domain,FDTD)方法是一种分析复杂介质电磁问题的得力方法.目前,国内外对于复杂介质电磁特性的分析主要集中在各向异性介质[1]、色散介质包括手征媒质[2]和铁氧体介质[3]及等离子体介质[4]等.但是对于双各向异性色散介质的研究,国内还没有看到报道,国外也还处于摸索起步阶段.双各向异性色散介质大多是新型的人工合成材料,可以应用到许多领域,如后向波媒体[5]、波导系统[6]以及微带天线的吸收器[7].由于其本构参数是张量的形式,导致分析这种介质变得相当复杂.为了研究该介质的电磁特性,需要对传统的FDTD方程进行拓展及改进.

基于分析色散介质电磁问题的改进FDTD方法,从现有文献看主要有辅助微分方程(Auxiliary Differential Equation, ADE)[8]、分段线性递归卷积方法(Piecewise Linear Recursive Convolution, PLRC)[9]、Z变换方法[10]和电流密度拉普拉斯变换时域有限差分(Current Density Laplace Transfer Finite-Difference Time-Domain, CDLT-FDTD)[11]方法.其中，Z变换方法在共振频率附近具有最高精确度,且容易编程实现.

基于一种改进的Z变换-FDTD方法,成功地研究了更一般情况下的双各向异性色散介质.Omega媒质[12]作为一种人工合成的双各向异性色散介质典型例子,便于建模分析,且具有一般代表性,因此选定该介质为基本研究对象.

1 双各向异性色散介质基本理论

如图1所示,形如Ω的金属粒子按部分平行xz和部分平行yz平面排列,称之为单轴型双各向异性色散介质.这些Omega媒质的共振行为导致了频率色散,使介质的本构参数变为频率的张量形式.

在频域下的双各向异性色散介质的本构方程可以写为：

(1)

(2)

式中：E、D、H、B分别是电场强度、电通量密度、磁场强度和磁通量密度；ε、μ、α是3×3阶的张量,分别称为张量介电系数、张量磁导率和张量电磁耦合系数,且都是ω的函数,它们的元素代表介质在各个方向上的本构参数； T表示转置.

(a) 主光轴沿波的传播方向

(b) 主光轴不沿波的传播方向图1 单轴型Omega介质板

1.1 单轴型双各向异性色散介质

根据图1(a)主光轴沿着z方向,单轴型双各向异性色散介质的本构参数为:

(3)

(4)

α(ω)=jK(ω)J.

(5)

假设在z方向Ω形金属粒子准静态极化,相对垂直的介电系数εn和磁导率μn可以忽略不计.则相对横向的本构参数分别可以表示为;

(6)

(7)

(8)

式中：ω0、ξ0、τK分别为谐振频率、阻尼系数和电磁耦合系数;εs、μs分别为静态相对介电系数和磁导率;ε∞、μ∞分别为频率远大于ω0时的相对介电系数和磁导率.

1.2 旋转形式

如图1(b)所示,主光轴与z轴有个夹角,设坐标轴以欧拉角α,β,γ旋转,则关于旋转坐标的双各向异性色散介质的本构参数张量形式变为:

εrotated(ω)=UT·ε(ω)·U;

(9)

μrotated(ω)=UT·μ(ω)·U;

(10)

αrotated(ω)=UT·α(ω)·U.

(11)

式中,U是旋转矩阵

(12)

2 双各向异性色散介质的FDTD算法

2.1 双各向异性色散介质本构方程的Z变换处理

由式(1)、(2)可知,Z域中的本构关系如下

D(z)=ε0ε(z)·E(z)T+

(13)

E(z)T.

(14)

式中,T是Z变换中离散时的采样周期.

将式(1)、(2)写成两部分,分别是与频率相关部分和频率无关部分,即

D(ω)=ε0(εc+εω(ω))·E(ω)+

(15)

B(ω)=μ0(μc+μω(ω))·H(ω)-

(16)

式中:

εc=ε∞Me+εnNe;

(17)

(18)

μc=μsMm+μnNm;

(19)

(20)

α(ω)=jK(ω)O.

(21)

对于旋转单轴型,有

(22)

其次,为了从频域转化到Z域,可以利用以下几个Z变换公式:

(23)

(24)

(25)

(26)

式中,β和ρ都是常量.则Z域中的本构参数张量形式ε(z)、μ(z)、α(z)可以写为:

(27)

(28)

(29)

2.2 双各向异性色散介质的修正场量方程

将ε(z)、μ(z)、α(z)代入Z域下的介质本构方程,得：

H(z);

(30)

(31)

式中:

(32)

(33)

(34)

注意,为了写法的简单和编程的方便,将式(30)、(31)写成以下形式:

D(z)=ε0εc·E(z)+ε0Se(z)z-1+

Sαh(z)z-1-Sαh(z)z-2;

(35)

B(z)=μ0μc·H(z)+μ0Sh(z)z-1-

2μ0Sh(z)z-2+μ0Sh(z)z-3-

Sαe(z)z-1+Sαe(z)z-2.

(36)

式中:

(37)

(38)

(39)

(40)

据此,对Z域下的修正场量方程进行移项化简,就可以得到修正的E(z)、H(z)方程:

Sαh(z)z-1+Sαh(z)z-2];

(41)

2μ0Sh(z)z-2-μ0Sh(z)z-3+

Sαe(z)z-1-Sαe(z)z-2].

(42)

已知Z域中有如下性质:z-nF(z)⟺f(t-nT).设定采样周期T等于FDTD算法中半个时间步长,则式(41)、(42)可写为:

(43)

(44)

在普通单轴型双各向异性色散介质情形下,有

(45)

(46)

结合(43)、(44)就可以推导出一维情形下修正的E、H方程.

2.3 双各向异性色散介质的修正S方程

将式(37)至(40)移项,并利用Z变换性质z-nF(z)⟺f(t-nT),可得：

(47)

(48)

(49)

(50)

3 数值验证及分析

3.1 一维单轴型双各向异性色散介质

计算厚度为10 cm的一层Omega板,所取参数为:ε∞=2,εs=5,εn=1.5,μ∞=1.1,μs=1.8,μn=1,f0=ω0/2π=2 GHz,ξ0=0.5,τK=0.7/ω0. 激励源选择高斯脉冲,计算结果如图2和图3所示,其中图2是反射波和透射波与解析解的对比图,图3是反射系数和透射系数与解析解的对比图.从两图可以看出,计算结果与解析解[13]相一致,证明了本方法的准确性.

3.2 一维旋转单轴型双各向异性色散介质

旋转单轴型的算法推导与第2节类似,仅在处理本构参数时引入旋转矩阵U即可.计算模型及参数与3.1节一样,欧拉角α=10°,β=15°,γ=45°.由于旋转矩阵U的影响,电磁耦合系数α变为一个复杂稠密的矩阵,导致了一个相对较强的电磁耦合性.因此,旋转单轴型的反射波和透射波就同时出现了同极化和交叉极化现象.图4是同极化下的反射波和透射波.图5是交叉极化下的反射波和透射波,它与同极化下的波形有较大区别,表现出了旋转矩阵的影响.另外,从两图可以看出计算结果与解析解[13]的误差很小,证明了Z-FDTD方法的可行性.

(a) 反射波

(b) 透射波图2 单轴型双各向异性介质的反射波和透射波波形

(a) 反射系数

(b) 透射系数图3 单轴型双各向异性介质的反射和透射系数

(a) 反射波

(b) 透射波图4 旋转单轴型双各向异性介质同极化反射波和透射波波形

(a) 反射波

(b) 透射波图5 旋转单轴型双各向异性介质交叉极化反射波和透射波波形

4 结 论

提出了一种适合于分析双各向异性色散介质电磁特性的方法,分别处理了普通单轴型和旋转单轴型两种不同情形的双各向异性色散介质的反射和透射情况,与解析解相符,表明了基于Z变换的修正FDTD方法的正确性.本文结果为将来解决三维实际双各向异性色散介质目标的电磁散射问题奠定了理论基础.

[1] 杨利霞, 葛德彪, 郑奎松, 等. 电各向异性介质FDTD并行算法的研究[J]. 电波科学学报, 2006, 21(1):43-48.

YANG Lixia, GE Debiao, ZHENG Kuisong, et al. Study of parallel FDTD algorithm for anisotropic medium on a PC cluster system[J]. Chinese Journal of Radio Science, 2006, 21(1):43-48.(in Chinese)

[2] PEREDA J A, GRANDE A, GONZALEZ O, et al. FDTD modeling of chiral media by using the Mobius transformation technique[J]. IEEE Antennas Wireless Propag Lett, 2006, 5(1):327-330.

[3] 胡振燕, 朱柏承, 周乐柱. 新型可调铁氧体-金属丝棋盘结构左手介质[J]. 电波科学学报, 2009, 24(5):799-803.

HU Zhenyan, ZHU Bocheng, ZHOU Lezhu. Novel tunable ferrite-wire chessboard structure left-handed metamaterial[J]. Chinese Journal of Radio Science, 2009, 24(5):799-803.(in Chinese)

[4] 杨利霞, 于萍萍, 马 辉, 等. 瞬变等离子体中电磁波频率漂移特性研究[J]. 电波科学学报, 2012, 27(1):18-23.

YANG Lixia, YU Pingping, MA Hui, et al. Frequency drifts characteristics for electromagnetic wave in suddenly creation plasma[J]. Chinese Journal of Radio Science, 2012, 27(1):18-23.(in Chinese)

[5] TRETYAKOV S A, SIMOVSKI C R, HUDLICKA M. Bianisotropic route to the realization and matching of backward-wave metamaterial slabs[J]. Phys Rev B, 2007, 75(15):3104-3109.

[6] YANG R, XIE J Y, LI X F, at el. Slow wave propagation in nonradiative dielectric waveguides with bianisotropic split ring resonator metamaterials[J]. Infrared Phys Technol, 2008, 51(6):555-558.doi:10.1016/j.infrared.2008.07.001.

[7] LANDY N I, BINGHAM C M, TYLER T, et al. Design, theory, and measurement of a polarization-insensitive absorber for terahertz imaging[J]. Phys Rev B, 2009, 79(12):5104-5109.

[8] NICKISCH L J, FRANKE P M. Finite-difference time-domain solution of Maxwell’s equations for the dispersive ionosphere[J]. IEEE Trans Antennas Propagat, 1992, 34(1):33-39.

[9] KELLEY D F, LUEBBERS R J. Piecewise linear recursive convolution for dispersive media using FDTD[J]. IEEE Trans Antennas Propagat, 1996, 44(6):792-797.

[10] YOUNG J L, NELSON R O. A summary and systematic analysis of FDTD algorithms for linearly dispersive media[J]. IEEE Antennas Propag Mag, 2001,43(1):61-126.

[11] 杨利霞, 谢应涛, 王祎君, 等. 一种适于1维磁等离子体电磁波传输特性的FDTD分析[J]. 强激光与离子束, 2009, 21(11):1710-1714.

YANG Lixia, XIE Yingtao, WANG Yijun, et al. Novel finite-difference time-domain analysis of electromagnetic wave transmission characteristics of magnetized plasma[J]. High Power Laser and Particle Beams, 2009, 21(11):1710-1714. (in Chinese)

[12] SAADOUN M M I, ENGHETA N. A reciprocal phase shifter using novel pseudochiral or Ω medium[J]. Microw Opt Technol Lett,1992,5(4):184-188.doi:10.1002/mop.4650050412.

[13] RIKTE S, KRISTENSSON G, ANDERSSON M. Propagation in bianisotropic media-reflection and transmission[J]. IEEE Proc Microw Antennas Propag, 2001,148(1):29-36.