牟春晖 陈娟2)† 范凯航 鲁艺

1) (西安交通大学信息与通信工程学院,西安 710049)

2) (西安交通大学深圳研究院,深圳 518057)

混合显-隐式时域有限差分方法的时间步长只由空间两个方向上的网格大小决定,因而该方法被广泛应用于沿一个方向具有精细结构的电磁目标的模拟.本文通过对传统时域有限差分方法基本公式进行近似,给出适用于线性、非色散空间,一维精细结构电磁目标模拟的通用混合显-隐式时域有限差分方法.在通用混合显-隐式时域有限差分方法中,通过媒质编号索引和迭代系数空间定位实现三对角矩阵的自动构建,引入等效参数方法解决麦克斯韦方程微分形式在不同介质分界面失效的问题,引入卷积完全匹配层边界条件对计算区域进行截断.此外,本文提出了一种缩减三对角矩阵数量的方法,该方法可有效降低程序运行所需内存,提高程序运行效率.应用通用混合显-隐式时域有限差分方法仿真了平面波照射介质板和双频微带倒F 天线两个模型,数值计算结果与传统时域有限差分方法和CST 软件的计算结果一致,但与传统时域有限差分方法相比,通用混合显-隐式时域有限差分方法的计算效率大大提高.

1 引言

时域有限差分(finite-difference time-domain,FDTD)方法[1,2]自提出以来已被广泛应用于电磁学研究的众多领域[3-9].传统FDTD 方法因采用显式差分算法而必须满足Courant-Friedrich-Levy时间稳定性条件[2].因此,传统FDTD 方法的时间步长受制于计算空间网格的最小值.当模拟目标在一个方向上具有精细结构时,为了保证计算精度,在目标精细结构所在的方向需要采用精细网格,依据时间稳定性条件,模拟此类问题时,时间步长需要取得很小,这就使得总的计算步数剧增.混合显-隐式FDTD(hybrid implicit-explicit FDTD,HIEFDTD)方法[10-12]通过对精细网格所在方向的空间求导项采用混合显-隐式差分算法,从而使得时间步长不再受限于精细网格大小.HIE-FDTD 方法自提出以来已被广泛应用于模拟在一个方向具有精细结构的问题,如基于石墨烯的各种电磁器件[13-16]、微带电路和天线[17,18]、线性/非线性集总网络[19,20]等.

在使用FDTD 方法计算空间任意一个网格上的电场和磁场时,需要用到该网格处介质的相关参数.因此在FDTD 方法的基本公式中,所有的迭代系数均需要标记空间位置.此外,由于微分形式的麦克斯韦方程在介质参数突变面失效,因此在编写FDTD 程序时需要引入等效参数[21,22]概念.然而,现有的关于HIE-FDTD 方法的文献和专著中所给出的HIE-FDTD 基本公式,大多数没有对迭代系数进行空间位置标记.依据这些公式编写的HIE-FDTD 程序,在仿真不同的电磁目标时,需根据目标介质分布情况手动构建三对角矩阵,这就导致HIE-FDTD 程序不具备通用性.少数文献在HIE-FDTD 迭代系数中增加了空间位置标记,但没有考虑介质分界面处微分形式麦克斯韦方程失效问题,在计算过程中也没有设置吸收边界,因此在模拟复杂目标时,计算精度难以保证,且不能计算开放区域问题[19,20].

本文提出了一种适用于线性、非色散空间,一维精细结构电磁目标模拟的通用HIE-FDTD 程序实现方法.在HIE-FDTD 基本公式中,对所有迭代系数进行空间位置标记,并对模拟空间中的所有媒质进行编号.在介质分界面引入等效参数概念.同时,引入卷积完全匹配层(convolutional perfectly matched layer,CPML)边界条件来截断计算区域.依据通用HIE-FDTD 方法编写的程序可根据目标介质分布情况自动构建三对角矩阵,且既可以模拟封闭区域问题,又可以模拟开放区域问题.此外,提出了一种缩减三对角矩阵数量的方法,该方法可有效降低程序运行所需内存,提高程序运行效率.本文应用通用HIE-FDTD 程序分别计算了平面波照射介质板时的透射场和反射场以及双频微带倒F天线S参数和辐射方向图,仿真结果与传统FDTD仿真结果以及CST 仿真结果相符合,且通用HIEFDTD 方法的计算效率远高于传统FDTD 方法.数值计算结果表明了该程序的通用性和有效性.

2 理论与程序实现

2.1 通用HIE-FDTD 基本公式

其中,C1,C2,C3和C4为FDTD 中电场和磁场时间推进公式中的迭代系数.若令

则(1)式可以重写为

假设目标精细结构位于z方向,则对空间求导项 ∂/∂z采用混合显-隐式差分.将(2)式进行如下变形:

其中,E1=(Az1-Az2+Bz1-Bz2).忽略(3)式的末尾项,得到

采用与传统FDTD 方法相同的Yee 元胞,并引入等效参数概念,将(5)式中的空间求导项采用二阶中心差分近似,便可得到电场和磁场各个场分量的迭代公式.为避免公式表述过于冗长,这里只给出了和两个场分量的迭代公式,其余4 个场分量的迭代公式可以参照这两个公式给出.

构建普惠金融体系成为当前我国新一轮农村金融改革的目标,小额信贷发展是实现普惠金融的重要内容。我国小额信贷发展过程中存在的问题主要体现在对“小额信贷”的理解存在偏差,小额信用贷款公司发展速度快，不仅存在“目标偏移”的倾向，而且社会资金的商业利益与社会责任脱节的情况也时有发生。指出我国小额贷款公司发展的总体目标是实现财务绩效与社会绩效协调发展,为实现我国普惠金融,在路径选择上应该注重农村金融市场机制培育。提出从开展社会绩效评价、降低交易成本等方面实现普惠金融的对策建议。

式中,Δx,Δy,Δz分别为沿x,y,z方向的空间网格大小;i,j,k为空间位置标记;n为时刻标记;Δt为时间步长.(6)式中的迭代系数分别为

其中,εx平均为等效介电系数,σx平均为等效电导率,µy平均为等效磁导系数,σmy平均为等效磁导率.等效参数的计算方法与传统FDTD 等效参数的计算方法相同[21,22].

使用FDTD 方法模拟开域电磁问题时,需要在计算区域的截断边界引入吸收边界条件.同理,使用HIE-FDTD 方法模拟开域问题,也需要引入吸收边界.由于采用CPML 作为吸收边界,电场和磁场分量均不需要分裂,计算公式较为简单,因此在通用HIE-FDTD 方法中仍采用CPML 边界条件作为吸收边界条件.在考虑CPML 吸收边界后,的计算公式变为

其中,

(8)式和(9)式中,Ψ为HIE-FDTD 中电场和磁场时间推进公式中卷积项的离散循环卷积形式,且CPML吸收边界各参数的取值与其在传统FDTD 方法中的取值一致.在非CPML 区域,kv的取值恒为1,av和bv的取值恒为0(v=x,y,z),此时(8)式退化为(6)式.引入CPML 后,电磁场其余场分量的迭代公式可参照(8)式给出[23-25].

显然,如(8)式所示,在计算n+1 时刻的Ex时需要用到同一时刻的Hy分量.将(8b)式代入(8a)式中,整理得到

其中,

同理,在计算n+1 时刻的Ey时需要用到同一时刻的Hx分量.采用与Ex相同的处理方法,得到Ey的计算公式如下:

其中,

综上所述,通用HIE-FDTD 方法的求解可通过如下方式完成: 首先,通过显式迭代方程计算;再利用(10)式和(11)式,通过解三对角矩阵方程求解;在此基础上,再通过显式迭代方程求解.完成一次时间步迭代,通用HIE-FDTD 方法共需要求解4 个显式迭代方程和2 个三对角矩阵方程.

2.2 程序实现

HIE-FDTD 方法与传统FDTD 方法在程序实现时最主要的差异在于,HIE-FDTD 需要构建三对角矩阵并求解其逆矩阵.以(10)式为例,在求解时,需要沿z方向以列为单位对其进行求解,且求解每一列时均需要构建一个三对角矩阵.假设计算区域沿x方向的网格数为Nx,沿y方向的网格数为Ny,沿z方向的网格数为Nz,且Nz=7,则求解每一列所需要构建的三对角矩阵方程为

其中,τ1k和τ2k分别为k取不同值时所对应的τ1和τ2的值,Exk为沿z方向第k个待求电场量,rk为k取不同值时(10)式等号右侧的已知量.显然,完成所有求解需要构建Nx×Ny个三对角矩阵,每个三对角矩阵的维度为(Nz—1)×(Nz—1),并执行Nx×Ny次的矩阵求逆计算.当计算区域较大时,执行程序将会耗费大量的内存和时间.为解决这一问题,本文提出了一种缩减三对角矩阵数量的方法.

以图1 所示模型为例,对通用HIE-FDTD 方法中三对角矩阵构建以及三对角矩阵数量缩减方法进行说明.如图1 所示,模型由目标和空气两部分组成,其中,目标由耶路撒冷十字金属层和长方体介质层组成.模型沿x和y方向的网格数均为18,沿z方向的网格数为10.由于本节只介绍三对角矩阵构建以及三对角矩阵数量缩减方法,因此这里只给出模型的网格尺寸,且整个模型用均匀网格线进行剖分,空间网格线如图1 中虚线所示.

图1 模型结构及网格分布Fig.1.Structure and grid distribution of the model.

对模型中的媒质进行编号,将空气标记为1,金属标记为2,介质标记为3.图2 为模型在x=4这一截面的介质分布矩阵.

图2 模型在x= 4 截面的介质分布矩阵Fig.2.Media distribution matrix of x= 4.

由图2 可见,当计算x=4 这一截面的时,共需要构建18 个三对角矩阵,每个三对角矩阵的维度为9×9.

当HIE-FDTD 方法中不引入等效参数概念时,三对角矩阵的值可由沿z方向每一列的介质组合决定.相同的介质组合对应的三对角矩阵是相同的.因此可以找出模型中不同的介质组合,求出这些不同介质组合对应的三对角矩阵即可.如图2 所示,虽然在x=4 这一截面有18 组介质组合,但多数介质组合是重复的,不重复的介质组合只有3 种,如图2 中阴影部分所示.

当HIE-FDTD 方法中引入等效参数概念时,参照传统FDTD 等效参数计算方法,D1,D2中的εx平均和σx平均是由(i,j,k),(i,j-1,k),(i,j,k-1)和(i,j-1,k-1) 4 个网格的ε和σ分别求平均得到,而D3,D4中的µy平均和σmy平均则是由(i,j,k)和(i,j-1,k)两个网格的µ和σm分别求平均得到.因此,对于相同的介质组合,其对应的三对角矩阵也不一定相同.如图2 中的y=3 和y=4两列的介质组合相同,但其对应的三对角矩阵的值却不同.因此,当在HIE-FDTD 方法中引入等效参数概念时,不可以再按照介质组合是否重复来缩减三对角矩阵数量.

由于三对角矩阵是基于τ1和τ2构建的,因此可以从τ1和τ2入手寻求减少三对角矩阵数量的方法.观察τ1和τ2的计算公式可以发现,τ1和τ2的值是随空间位置变化的.因此可以先计算出空间所有网格点的τ1和τ2,然后将其按列(沿z方向)分别存储到二维矩阵T1和T2中.在此基础上,令T=,则T矩阵中的每一列对应求解时所需要构建的一个三对角矩阵.T1矩阵、T2矩阵和T矩阵的结构如图3 所示,其中,τ1-i,j,k表示空间位置(i,j,k)处的τ1值,τ2-i,j,k表示空间位置为(i,j,k)处的τ2值,Nz为计算区域沿z方向的网格数.

图3 T1 矩阵、T2 矩阵和T 矩阵结构图Fig.3.Structure of T1,T2 and T.

根据电磁目标的结构特征,可以预见,T矩阵中可能存在一些相同的列,基于这些相同的列所构建的三对角矩阵是完全相同的.因此可以在构建三对角矩阵之前,先剔除T矩阵中的重复列,这样便可以有效减少三对角矩阵的数量.在每一个时间步中计算时,只需要根据空间位置标记索引到对应的三对角矩阵进行求解即可.

图1 所示模型结构简单,网格规模较小.实际中的模型往往结构复杂,尺寸较大[26].为说明三对角矩阵缩减方法的有效性,以图1 模型中的目标部分为基础,分别构建了5×5 和10×10 的目标阵列,并在目标阵列外围各设置2 层空气.5×5 目标阵列的网格规模为74×74×10,10×10 目标阵列的网格规模为144×144×10.分别采用直接法和缩减法对单目标模型和阵列目标模型进行三对角矩阵构建及其逆矩阵计算.两种方法需要构建的三对角矩阵数量、程序运行时间以及程序所需内存如表1 所列.由表1 可见,采用缩减法可以有效减少三对角矩阵数量,有效降低程序运行所需内存,有效缩短程序运行时间.

表1 直接法和缩减法对比结果Table 1.Comparison results of the direct method and the reduction method.

3 数值结果

为验证本文所提出的通用HIE-FDTD 程序实现方法的有效性,依据该方法编写了Matlab 程序,并用程序计算了两个数值算例,分别为平面波照射介质板和双频微带倒F 天线.在这两个例子中,目标的精细结构均位于z方向,因此对于空间求导项∂/∂z采用混合显-隐式差分算法.

3.1 平面波照射介质板模型

介质板由介质基板和金属表面涂层两部分组成.介质基板的介质参数为εr=2.2,µr=1,σ=0.04,σm=0 .介质基板沿x和y方向的长度均为210 mm,沿z方向的厚度为0.2 mm.金属涂层由5×5 个耶路撒冷十字组成,厚度为0.2 mm.介质板结构及耶路撒冷十字尺寸如图4 所示.整个计算空间网格数为120×120×56.网格大小分别为:在空气中,Δx=Δy=Δz=5 mm;在介质基板和金属涂层中,Δx=Δy=3 mm,Δz=0.1 mm.

图4 介质板结构及尺寸Fig.4.Structure and dimensions of the dielectric plate.

假设平面波从金属涂层一侧照射到介质板上,在介质板中轴线上设置两个观察点,计算观察点处电场的时域波形.观察点1 位于介质板有金属涂层一侧,且距金属涂层的距离为10.2 mm.观察点2位于介质板无金属涂层一侧,且距无金属涂层的距离为10.2 mm.平面波的激励源为高斯脉冲,其时域表达式为

其中,τ=3.5×10-10,t0=0.8τ.

观察点1 处Ex的时域波形如图5(a)所示,观察点2 处Ex的时域波形如图5(b)所示.从图5 可以看到,通用HIE-FDTD 程序的计算结果与传统FDTD 程序的计算结果相符合.从仿真时间上来看,传统FDTD 方法的仿真时间为300 min,而通用HIE-FDTD 方法的仿真时间仅需32 min,约为传统FDTD 方法仿真时间的1/10.此外,采用直接法的HIE-FDTD 程序的仿真时间为45 min,与通用HIE-FDTD 程序仿真时间的差异主要是因为需构建更多三对角矩阵及其逆矩阵造成的.可以预见,当模型网格规模更大时,采用直接法的HIEFDTD 程序的计算时间将会进一步增加.

图5 观察点处Ex 的时域波形 (a) 观察点1;(b) 观察点2Fig.5.Time domain waveforms of Ex at the observation point: (a) Observation point one;(b) observation point two.

3.2 双频微带倒F 天线模型

双频微带倒F 天线的结构如图6 所示.金属地板和倒F 结构分别刻蚀在介质板的两侧.介质基板的介质参数为εr=2.2,µr=1,σ=0,σm=0.介质基板沿x和y方向的长度分别为38.4 和39 mm,沿z方向的厚度为0.8 mm.整个计算空间网格数为56×54×46,且沿x和y方向的最小网格值为2.4 mm,沿z方向的最小网格值为0.2 mm.

图6 双频微带倒F 天线结构图Fig.6.Structure of the dual-frequency microstrip inverted-F antenna.

双频微带倒F 天线的S11仿真结果如图7 所示.可以看出,倒F 天线的工作频率为2.43 和4.97 GHz.图7 同时给出了传统FDTD,HIE-FDTD 和CST仿真结果,三条S11曲线符合度较好.图8 给出了天线在2.43 GHz 的辐射方向图,图9 给出了天线在4.96 GHz 的辐射方向图.从图8 和图9 可以看出HIE-FDTD 的仿真结果与CST 仿真结果基本一致.从仿真时间上来看,传统FDTD 方法的仿真时间为420 s,通用HIE-FDTD 方法的仿真时间为90 s,约为传统FDTD 方法仿真时间的1/5.此外,采用直接法的HIE-FDTD 程序的仿真时间为113 s.可见当模型网格规模较小时,缩减三对角矩阵数量方法对于总仿真时间的减少作用并不明显.此时,仿真时间主要由时间步长和仿真步数决定.

图7 倒F 天线S11Fig.7.S11 of the inverted-F antenna.

图8 天线在2.43 GHz 的辐射方向图 (a) XOY 面;(b) XOZ 面Fig.8.Radiation patterns of the inverted-F antenna at 2.43 GHz: (a) XOY plane;(b) XOZ plane.

图9 天线在4.96 GHz 的辐射方向图 (a) XOY 面;(b) XOZ 面Fig.9.Radiation patterns of the inverted-F antenna at 4.96 GHz: (a) XOY plane;(b) XOZ plane.

4 总结

本文提出了一种通用HIE-FDTD 方法,并依据该方法编写了一套通用HIE-FDTD 程序.使用该套程序可以实现对线性、非色散空间中,在一维方向具有精细结构的任意电磁目标的仿真.通用HIE-FDTD 程序可实现三对角矩阵的自动构建,并对三对角矩阵数量进行自动缩减,使用一套程序即可实现对任意目标的模拟,而不是像现有HIEFDTD 程序那样仿真一个模型便需要修改一次程序.另外,在通用HIE-FDTD 方法中,计算迭代系数时考虑了微分形式麦克斯韦方程在介质分界面失效的问题,并引入等效参数方法解决这一问题.在以往的HIE-FDTD 文献和专著中均没有提到等效参数方法这一概念.在通用HIE-FDTD 基本公式中,加入CPML 吸收边界项,使得通用HIEFDTD 程序不仅可以模拟封闭区域问题,也能模拟开放区域问题.使用通用HIE-FDTD 程序仿真了平面波照射介质板模型以及双频微带倒F 天线模型,仿真结果与传统FDTD 仿真结果以及CST软件仿真结果相符合,且通用HIE-FDTD 程序的仿真时间远少于传统FDTD 的仿真时间.本文提出的通用HIE-FDTD 方法只适用于非色散空间,后续研究中可结合辅助方程等方法将其扩展到色散空间,使其应用范围更加广泛.