Padé近似下模拟一般色散媒质的FDTD改进方案

刘广东

·物理电子学·

Padé近似下模拟一般色散媒质的FDTD改进方案

刘广东

(阜阳师范学院物理与电子工程学院 安微 阜阳 236037)

为了方便模拟不同媒质的色散特性，提出了一种时域有限差分(FDTD)改进方案，适用于统一处理几类各向同性、线性、有磁电色散媒质的电波传播问题：媒质类型可以是Havriliak-Negami(H-N)媒质、Davidson-Cole(D-C)媒质、Cole-Cole(C-C)媒质、Debye媒质、常规(非色散)媒质或其任意组合；媒质属性可以是单极或多极的、有电耗的或无电耗的。该方案利用帕德(Padé)近似法，导出了一组整数阶的辅助微分方程(ADEs)，既克服了其中分数阶导数的主要困难，又展现了通用性好、复杂度低的优势。通过对一维及三维算例解析、数值结果之间的对比，初步证实了改进方案的可行性和有效性。

辅助微分方程; 时域有限差分法; 色散媒质; 帕德近似法

实验测量证实：生物组织、大地、岩石、水、液晶、聚合物材料等众多媒质的电参数展现了频率相关性(这些媒质称作电色散媒质)；不同媒质还可能展现不同的电色散特性[1-2]。至今已提出好几种经验模型以描述不同的电色散特性，一般可概括为两大类型[3]：第一大类主要有Debye模型[4](如人体组织、等离子体、土壤)、Drude模型[5](如等离子体、金属)和Lorentz模型[2](如光学材料、人工介质)，其共性在于：媒质的(复)介电常数是jω(j表示虚数单位，ω表示角频率)整数次幂的函数[3]；第二大类主要包含Cole-Cole(C-C)模型[6-8](如生物组织、高分子材料)、Davidson-Cole(D-C)模型[3](如甘油、丙二醇)和Havriliak-Negami(H-N)模型[9](如聚合物)，其共性在于：媒质的(复)介电常数是jω分数次幂的函数[3]。后文将以不同媒质对应的色散模型(或所属大类)命名以区分不同的媒质类型(如称色散特性满足C-C模型的媒质为C-C媒质)。

近几年来，由于媒质的色散特性引起了较多关注[10-11]，相应地，FDTD法的解算对象已由常规的非色散媒质[12]拓展到色散媒质[6,13]。

FDTD建模第一大类(Debye、Drude、Lorentz)媒质，已发展了多种方法，主要有3种思路：1) 引入ADE[2,7-9]；2) 利用递归卷积(recursive convolution, RC)[2,13]；3) 定义移位算子(shift operator, SO)[2,14]。

然而，FDTD建模第二大类(C-C、D-C、H-N)媒质，仍是当前的难点之一，其主要困难是差分实现分数阶导数[3]。对此，文献[15]进行了一系列开创性的基础工作：利用帕德(Padé)近似法，导出了一组整数阶ADEs，克服了分数阶导数困难，分别提出了处理(无磁、线性、各向同性的)单极、多极C-C媒质[7-8]、单极D-C媒质[3]和单极H-N媒质[9]的FDTD方案；一维的系列仿真结果初步证实了这些方案的可行性[7-9]。文献[6]给出了建模多极C-C媒质的FDTD改进方案和三维的仿真测试。

为了增加通用性和降低复杂度，方便统一编程处理不同类型色散媒质的电波传播问题，本文主要以文献[7-9]为基础，提出了一种FDTD改进方案，创新(改进)之处主要体现在：1) 从单极色散媒质推广到多极情形；2) 从无磁媒质推广到有磁情形；3) 增加了静态电导率项；4) 增加了三维数值算例检验该方案的性能和精度。

1 问题描述

首先给出后文满足的假设条件和必要的符号说明：1) 假设所研究的色散媒质是各向同性、线性、有磁和电色散的，该假设条件经常应用于地球物理勘探、微波医学成像、红外与毫米波技术等领域[1-5]；2) 符号顶部加点表示其值为复数。

若问题空间r处存在一般的(同时存在电极化损耗和欧姆损耗)多极H-N媒质，则频域本构关系为：

式中，D˙和E˙分别表示电位移矢量和电场强度矢量；ω = 2πf，f为频率；ε0为真空介电常数；多极H-N媒质的相对介电常数定义为[16]：

式中，ε∞为光学相对介电常数；σs为静态电导率；Δεw= εs,w−ε∞；εs,w为第w极的静态相对介电常数；W为极总数；第w极的复值因子wf˙为：

式中，αw(0≤αw≤1)；βw(0≤βw≤1)为第w极色散参数；wz˙= jωτw(r)；τw表示第w极的弛豫时间。

本文采用的色散模型是较为通用的形式：1) 当 W = 1时，简化为单极情形；2) 当W = 1且αs= 0时，简化为文献[9]的式(1)；3) 如果W = 1、αw= 1且σs= 0时，退化为文献[3]的式(1)；4) 如果βw= 1且σs= 0时，退化为文献[8]的式(1)；5) 如果W = 1、βw= 1且σs= 0时，变为文献[7]的式(1)；6) 如果αw= 1且βw= 1时，变为Debye模型[4]；7) 假如αw= 1、βw= 1且Δεw= 0时，转化为非色散情形。

2 辅助微分方程(ADEs)

为了克服FDTD方案中分数(αw、βw一般为分数)阶导数的困难，首先，借鉴文献[7-9]的研究思路，利用帕德(Padé)近似法[15]，并由wf˙在wz˙= 1处解析，将wf˙近似为：

式中，N、M分别表示分子、分母多项式的阶数；pn、qm分别表示分子、分母多项式的系数，可通过求解线性方程组获得，相关细节参见文献[7-9]。

将式(3a)代入式(2a)，易得相对介电常数的帕德近似值为：

这样，相对介电常数已全部转化为jω整数次幂的函数，有利于后文FDTD方案的讨论。

其次，定义第w极(辅助的)频域电极化强度矢量为：

最后，应用傅里叶逆变换(inverse Fourier transform, IFT)，可以得到一组(W个)关于时域电极化强度矢量Pw的ADEs为：

其中：

式中，E表示时域电场强度矢量；t表示时间。至此，这组辅助微分方程由于包含了整数阶的时间导数，可以方便地通过有限差分法实现时间的离散化。

3 时域有限差分方案

前述电色散媒质中，时域(微分形式)的麦克斯韦(Maxwell)方程组为[8]：

式中，H(为了表达简洁，省略了自变量，下文同)表示磁场强度矢量H(r,t)；μr表示相对磁导率，当μr= 1时，媒质简化为文献[7-9]的无磁特例。

在t=iΔt(i和Δt分别表示时间步的指标和步长)时刻的电场强度矢量Ei，其n阶的时间偏导数可通过中心差分近似为[7-9]：

式中，l=n/2表示对n/2上取整；Qnm为中心差分参数[6-9]；Pw的中心差分近似可类似地获得。至此，就实现了式(5)的FDTD时间离散化。

文献[9]显示，欲使FDTD方案稳定，需满足N≤M。为了简化描述，假定M = N(其原因后文给出)，则有lM=lN。进一步假定已知iΔt时刻的Pw、E和(i+0.5)Δt时刻的H，分别差分式(8)～式(9)，并引入辅助矢量Φw[8]，经过整理，多极H-N媒质的FDTD时间迭代步进方案可简述如下：

1) 计算辅助矢量为：

其中：

式中，k = −lM,−lM+1,…,lM；ν = max{2|k|−1,0}。

2) 更新(i+1)Δt时刻的电场强度矢量为：

3) 更新(i+1)Δt时刻的电极化强度矢量为：

4) 更新(i+1.5)Δt时刻的磁场强度矢量为：

前述FDTD方案中，对于每个离散的网格，每个场分量的存储量为M(W+1)+W，磁场以外的计算时间(乘法运算次数)为(2M+1)W+3，当W=1时与文献[9]的结果相同。

该FDTD方案的时间迭代步进流程(技术路线)如图1所示[6]。

图1 FDTD时间迭代步进流程图

4 结果与讨论

由于前述FDTD方案中采用了帕德(Padé)近似，下面首先检验帕德近似的精度。

4.1 帕德近似的精度

定义帕德近似的相对均方(relative mean square, RMS)误差为[3,7]：

式中，媒质相对介电常数的解析值ε˙r和Padé近似值ε˙ˆr可分别通过式(2a)和式(3b)求得；fH、fL分别表示上、下限频率。

算例 1 以文献[9]中的均匀单极H-N媒质(聚丙烯酸正辛酯，聚丙烯酸的一种衍生物)为测试目标，各模型参数分别为W = 1、ε∞= 2.61、σs= 0.001 S/m、εs,1= 3.88、τ1= 39.98 ps、α1= 0.73、β1= 0.66。上、下限频率分别取为fH= 100 GHz、fL= 10 MHz，分母多项式的阶数取为M = 3,4,…,10，分子多项式阶数分别取为N = M、N = M−1和N = M−2时，帕德近似的相对均方误差e随M的变化关系如图2所示(误差轴为对数坐标，后文同)，其中当M=3、4、5、6时，本文的误差与文献[9]的误差对比如表1所示。

图2 算例1相对均方误差随M的变化关系

表1 相对均方误差随N、M的变化关系

分析图2发现：1) 当M相同时，在N = M、N = M−1、N = M−2的3种情形下，第一种情形误差最小，这正是前文讨论FDTD时间步进方案时选取N = M的依据；2) 在N = M、N = M−1、N = M−2的3种情形下，随着M的增加，误差均逐渐减小(然而，实现FDTD方案所需的存储量M(W+1)+W和计算时间(2M+1)W+3也同时增加，建议折中选择M、N)。

结合表1可知：1) 在表中所列举的12种情形下，本文的误差均略大于文献[9]的误差(适当增加离散频点数可使二者相当)，因此，文献[9]的单极情形可视为本文改进方案的一个特例；2) 当N = M = 4时，本文的误差约为0.003 6，不仅计算精度能被一般的工程应用所接受，而且计算速度比较适中。

算例 2 (为不失一般性)选取一种2极(即W=2)的均匀H-N媒质为测试目标，其余模型参数分别为ε∞= 2.61、σs= 0.01 S/m、εs,1= 3.88、εs,2= 3.50、τ1= 39.98 ps、τ2= 29.98 ns、α1= 0.73、α2= 0.80、β1= 0.66、β2= 0.70。频率上、下限及多项式的阶数均与算例1相同，相对均方误差e随M的变化关系如图3所示。

由图3可见：对于多极H-N媒质，也有类似于算例1(单极情形)的结果，只是误差略微增加。当N = M = 4时，相对均方误差约为0.004 3，还能满足一般的工程应用需求，因此，在后文的FDTD方案中，多项式的阶数也默认这组取值。4.2 FDTD方案的可行性

图3 算例2相对均方误差随M的变化关系

为了检验前述FDTD改进方案的可行性，再给出两个算例，选择算例基于以下考虑：1) 问题空间涵盖不同物理维度；2) 色散媒质涵盖不同类型；3) 媒质涵盖不同极数；4) 色散模型参数涵盖文献作者自拟和真实测量；5) 测试采用不同类型目标参数。

算例 3 出自文献[3]的一维问题：单极均匀无磁D-C媒质(媒质参数分别为W = 1、μr= 1、ε∞= 2、σs= 0 S/m、εs,1= 50、τ1= 153 ps、β1= 0.9和α1= 1.0)充满z≥0半空间，其余空间为真空；一列入射平面波(电场沿x方向极化)沿+z方向传播，平面波由下述的超宽带(ultra-wideband, UWB)调制高斯脉冲源s激励产生[3]：

式中，参数a =1.26×1010s−1，中心频率fc= 6 GHz，该脉冲的频谱主要涵盖0.1～10 GHz频率范围。

FDTD离散化的空间、时间步长分别取为Δz = 1.1 mm、Δt = 1.5 ps，离散网格数、时间步数分别取为nz= 200、nt= 3 000。仿真区域−z一侧采用简单的吸收边界：Ex(zmax, t+Δt) = Ex(zmax−Δz, t)，其中zmax表示对应的边界位置坐标。

记z、z+d分别表示D-C媒质中的两个不同位置，FDTD仿真并且存储这两处的时域电场Ex(z, t)、Ex(z+d, t)，后处理(傅里叶变换)获得相应的频域电场E˙x(z, ω)、E˙x(z+d, ω)，则色散媒质相对介电常数的FDTD估计值为[3,7]：

式中，c0表示真空光速；参数γR、γI分别为：

选取d = 20Δz时，复值相对介电常数实部和负虚部的解析值、Padé近似值、FDTD估计值随频率的变化关系如图4所示(频率轴为对数坐标，后文同)。

从图4可以看出，在UWB条件下，以一维问题的相对介电常数为测试目标，其FDTD估计值和Padé近似值或解析值之间差异不大(实际结果显示相对均方误差均小于1%)；再与文献[3]中图4a的结果对比发现，二者也基本保持一致，因此，文献[3]的单极情形可视为本文改进方案的又一个特例。

图4 复值相对介电常数随频率的变化关系

算例 4 选取一个更具代表性的三维问题：计算空气(视为真空)中均匀4极C-C弱磁性介质球的后向雷达散射截面(radar cross section, RCS)[2]，球半径为50 mm，由浸润型脂类组织构成，其模型参数(源于实验测量[17])分别为ε∞= 2.5、σs= 0.035 S/m、Δε1= 9.0、Δε2= 35、Δε3= 3.3×104、Δε4= 1.0×107、τ1= 7.96 ps、τ2= 15.92 ns、τ3= 159.15 μs、τ4= 15.915 ms、β1= 0.80、β2= 0.90、β3= 0.95、β4= 0.99；其余模型参数分别为μr= 1.2、W = 4、α1= α2= α3= α4= 1.0；采用的入射平面波与算例3相同。

FDTD离散化的空间网格尺寸、时间步长分别取为Δ = Δx= Δy= Δz= 3 mm、Δt = Δ/(2.c0)，离散获得的网格数、时间步数分别为n = 50×50×50、nt= 5 000。主FDTD区域周围用10层卷积完全匹配层(convolution perfectly matched layer, CPML)[2]吸收边界截断。

作为比较，RCS的解析值由Mie级数法[2]获得，两种方法的计算结果对比如图5所示。

通过图5发现，在较宽泛的频率范围，以三维吸波介质球的后向RCS为测试目标，其FDTD估计值和解析值仍然符合较好(结果显示，相对均方误差较一维情形有所增加，原因可能是网格剖分粗糙所致，有望通过先进的网格剖分技术进一步减小误差[2])。

综合以上4个算例可看出，利用前述FDTD改进方案模拟UWB电磁波在C-C、D-C、H-N等几类电色散媒质中传播，估计相对介电常数、后向雷达散射截面所产生的误差，一般的工程应用可以接受。这些仿真结果初步证实，利用本文发展的FDTD方案，模拟一般电色散媒质中的电波传播是可行的。

图5 介质球后向雷达散射截面随频率的变化关系

5 结 束 语

本文改进了文献[7-9]提出的FDTD方案，改进后的方案具有适用范围广、实现复杂度低等优势，适用于统一处理射频、微波、红外和毫米波等工程中一般电色散媒质的电波传播问题，媒质可以是一维或多维的、单极或多极的、无磁或有磁的。算例的初步结果显示，在超宽带频谱范围，该方案具有较高的数值精度。该方案的鲁棒性还需经历实验测量、噪声影响等因素的检验。

[1] 陈西良, 陈欣, 朱智勇. 不同极性聚合物材料的THz-TDS光谱测量研究[J]. 红外与毫米波学报, 2013, 32(2): 150-159. CHEN Xi-liang, CHEN Xin, ZHU Zhi-yong. THz-TDS spectra study of polymer materials with different polarity[J]. Journal of Infrared and Millimeter Waves, 2013, 32(2): 150-159.

[2] 葛德彪, 闫玉波. 电磁波时域有限差分法[M]. 第3版. 西安: 西安电子科技大学出版社, 2011: 259-300. GE De-biao, YAN Yu-bo. Finite-difference time-domainmethod for electromagnetic waves[M]. 3rd ed. Xi’an: Xidian University Press, 2011: 259-300.

[3] REKANOS I T. FDTD schemes for wave propagation in Davidson-Cole dispersive media using auxiliary differential equations[J]. IEEE Transactions on Antennas and Propagation, 2012, 60(3): 1467-1478.

[4] 刘广东, 张业荣. 一种处理分层有耗色散介质的时域逆散射方法[J]. 电子学报, 2011, 39(12): 2856-2862. LIU Guang-dong, ZHANG Ye-rong. An approach to the time-domain inverse scattering problem for the stratified frequency-dispersive lossy media[J]. Acta Electronica Sinica, 2011, 39(12): 2856-2862.

[5] WANG Ai-hua, CAI Jiu-ju. Modeling of radiative properties of metallic microscale rough surface[J]. Journal of Central South University, 2012, 19: 1482-1487.

[6] 刘广东, 张开银, 范士民. 一种处理Cole-Cole色散媒质的FDTD改进方案[J]. 计算物理, 2014, 31(2): 257-265. LIU Guang-dong, ZHANG Kai-yin, FAN Shi-min. A modified FDTD scheme for wave propagation in general Cole-Cole dispersive media[J]. Chinese Journal of Computational Physics, 2014, 31(2): 257-265.

[7] REKANOS I T, PAPADOPOULOS T G. An auxiliary differential equation method for FDTD modeling of wave propagation in Cole-Cole dispersive media[J]. IEEE Transactions on Antennas and Propagation, 2010, 58(11): 3666-3674.

[8] REKANOS I T, PAPADOPOULOS T G. FDTD modeling of wave propagation in Cole-Cole media with multiple relaxation times[J]. IEEE Antennas and Wireless Propagation Letters, 2010(9): 67-69.

[9] REKANOS I T. FDTD modeling of Havriliak-Negami media[J]. IEEE Microwave and Wireless Components Letters, 2012, 22(2): 49-51.

[10] 刘广东, 张业荣. 二维有耗色散介质的时域逆散射方法[J]. 物理学报, 2010, 59(10): 6969-6979. LIU Guang-dong, ZHANG Ye-rong. Time domain inverse scattering problem for two dimensional frequencydispersive lossy media[J]. Acta Physica Sinica, 2010, 59(10): 6969-6979.

[11] 曹苗苗, 刘文鑫, 王勇, 等. 介质加载复合光栅结构的色散特性研究[J]. 物理学报, 2014, 63(2): 024101. CAO Miao-miao, LIU Wen-xin, WANG Yong, et al. Dispersion characteristics of dielectric loaded metal grating[J]. Acta Physica Sinica, 2014, 63(2): 024101.

[12] 刘瑜, 梁正, 杨梓强. 混合并行技术在FDTD计算中的应用研究[J]. 电子科技大学学报, 2009, 38(2): 222-226. LIU Yu, LIANG Zheng, YANG Zi-qiang. Study and application on hybrid parallel FDTD algorithm[J]. Journal of University of Electronic Science and Technology of China, 2009, 38(2): 222-226.

[13] 王飞, 魏兵. 电各向异性色散介质电磁散射时域有限差分分析的半解析递推卷积方法[J]. 物理学报, 2013, 62(4): 044101. WANG Fei, WEI Bing. Semi-analytical recursive convolution algorithm in the finite-difference time domain analysis of anisotropic dispersive medium[J]. Acta Physica Sinica, 2013, 62(4): 044101.

[14] 魏兵, 董宇航, 王飞, 等. 基于移位算子时域有限差分的色散薄层节点修正算法[J]. 物理学报, 2010, 59(4): 2443-2450. WEI Bing, DONG Yu-hang, WANG Fei, et al. A modificatory algorithm for electrically thin dispersive layers base on shift operator finite-difference time-domain method[J]. Acta Physica Sinica, 2010, 59(4): 2443-2450.

[15] BAKER G A, GRAVES M P. Padé approximants[M]. New York, USA: Cambridge University Press, 1996: 1-8.

[16] HAVRILIAK S, NEGAMI S. A complex plane representation of dielectric and mechanical relaxation processes in some polymers[J]. Polymer, 1967, 8(4): 161-210.

[17] GABRIEL S, LAU R W, GABRIEL C. The dielectric properties of biological tissues III. Parametric models for the dielectric spectrum of tissues[J]. Physics in Medicine and Biology, 1996, 41: 2271-2293.

编 辑 黄 莘

Improved FDTD Scheme Based on Padé Approximant Method for Modeling of Wave Propagation in General Dispersive Media

LIU Guang-dong

(School of Physical and Electronic Engineering, Fuyang Normal College Fuyang Anhui 236037)

A modified finite-difference time-domain (FDTD) scheme is developed to simulate wave propagation in different electrically dispersive media with isotropic, linear and magnetic properties. The presented scheme is applicable to several types of general frequency-dependent media such as Havriliak-Negami (H-N), Davidson-Cole (D-C), Cole-Cole (C-C), Debye dispersive media or nondispersive media, which are lossless or lossy, with single pole or multiple relaxation times. The main difficulty in this scheme is the appearance of fractional derivatives. Based on the Padé approximant method, a set of auxiliary differential equations (ADEs) of integer order are derived. Thus, this difficulty is circumvented, and its advantage in universality and complexity is also exhibited. The feasibility and validity of the presented scheme are preliminarily demonstrated by the comparisons between analytic and numerical results from several one-dimensional (1-D) and three-dimensional (3-D) examples.

auxiliary differential equation (ADE); finite-difference time-domain (FDTD) method; frequency-dependent media; Padé approximant method

O441.4

A

10.3969/j.issn.1001-0548.2015.06.009

2014 − 05 − 30；

2015 − 06 − 11

国家自然科学基金(51271059)；安徽省教育厅自然科学研究重点项目(KJ2014A193)；安徽省科技计划(12010302080, 1501031114)

刘广东(1972 − )，男，博士，副教授，主要从事电磁散射和逆散射方面的研究.