二次B-Spline时域基函数的TDFEM的应用

吴 霞，周乐柱

(1.北京大学信息科学技术学院，北京 100871，wuxiaeecs00@gmail.com;2.复旦大学信息科学与工程学院，上海 200433)

二次B-Spline时域基函数的TDFEM的应用

吴 霞1，2，周乐柱1

(1.北京大学信息科学技术学院，北京 100871，wuxiaeecs00@gmail.com;2.复旦大学信息科学与工程学院，上海 200433)

提出了时域有限元法(TDFEM)的一种新的基函数—二次B-spline时域基函数.首先简述了时域有限元法的原理和基本公式;然后提出了新型的基于B-spline函数的条件稳定和无条件稳定的时域有限元法方案，并应用于三维电磁辐射问题.通过典型的算例对这两种方案的精度、运算时间进行了比较，证实了基于二次B-spline函数的时域有限元法的有效性.通过稳定性理论分析得出该算法的精确稳定性，并且通过数值计算的结果得到验证.

时域有限元法;二次B-spline时域基函数;电磁辐射;无条件稳定;超宽带

在计算电磁领域，有不同的时域分析计算方法，例如时域有限差分法(FDTD)、时域积分法(TDIE)和时域有限元法(TDFEM)等，致力于对各种不同的结构进行建模与仿真分析.在众多的方法中，TDFEM因其继承了FEM特别适用于复杂的几何结构和介电特性分布的优点，及在非均匀介质体内产生的也是稀疏矩阵的长处，因而在解决宽带特性的电磁场问题上弥补了FDTD和TDIE法的某些不足而被人们渐渐关注.TDFEM起源于上世纪80年代，分为两大类，一类是最初由 A.C.Cangellaris等［1－3］提出的微分形式的点匹配时域有限元法.该方法结合了FDTD显式积分法的简明性和传统有限元的灵活性，但是会出现源于电场和磁场点在网格上的交替出现而产生的类似FDTD法的越级(leap-frog)问题.为了避免或解决这些问题，第2类隐式时域有限元法在1995年提出［4］.该方法是基于电场(或磁场)的矢量有限元空间基函数和时间基函数展开的二阶矢量波方程，其优点是简单而稳定，缺点是每个时间步必须求解一个矩阵方程;对于开域空间问题，采用传统的吸收边界条件，有限元计算空间十分庞大，这个问题更加严重，以致无法进行.由于这些原因，这类时域有限元法发展缓慢.直到2000年初，由于完全匹配层(PML)［5］等概念的提出和应用，时域有限元法才得到较大的发展.

本文研究上述第2类的时域有限元法，与已报道的TDFEM方法不同的是，本文首次引入了B-spline函数作为时间基函数.另外，利用Newmark-β方法验证了相应的检测函数.B-spline函数早在1946年被I.J.Schenberg提出［6］，因其插值特性和计算简单的优点近期才在计算电磁算法中被重视.B-spline时域基函数已经在FDTD和TDIE算法中被采用并取得了不错的效果［7－8］.然而，目前在时域有限元法的研究中还没有具体应用B-spline时域基函数的报道.本文算法的另外两个特点是:1)采用Sacks PML做截断边界条件;2)采用高阶矢量有限元作为空间基函数.

1 TDFEM基本公式

为了推导出电场的时域有限元方法，从一个时域电磁场的初值问题开始［3］，方程式如下:

其中J0代表外加电流源，εr和μr分别是相对介电常数和相对磁导率，c0代表自由空间光速.用Sacks PML来截断有限区域，在PML区域，电场满足下式:

为了把式(4)的半离散化的时域方程转变成完全离散化的时域方程，ej(t)应当由某种时间基函数Tm(t)展开，即

其中emj是ej(t)在t=tm(m为整数)时相对应的值.在以往报道的时域有限元法文献中，Tm(t)都采用如下的线性时域基函数:

本文采用二次B－spline函数作为时域基函数.

2 B-spline时域基函数

构造B样条函数时，给定阶数n，B样条函数φn(x)的表达式为［11］

其中:H(x)是阶跃函数.B样条函数虽然是分段函数，但(n－2)次导数连续，因此用在TDFEM中可以更精确地进行时域离散.

在式(6)中取n=3得到本文所用的二阶B-spline时域基函数的表达式由下式给出.

3 两种基于B-spline的TDFEM

该新算法的稳定性是至关重要的.鉴于稳定性依赖于时域的离散化，而且有多种方法把式(4)变成一个完全离散的矩阵方程，如中心差分、前向差分、后向差分或Newmark-β 法［12］.本文采取和空间离散同样的方法，即用时域基函数来获得完全离散化的方程.时域离散展开系数的值取决于离散化时域基函数和测试函数的具体形式.

3.1 条件稳定(CS)B-spline函数TDFEM方案

第一种选择是采用如下的B-spline函数和测试函数:

将时域基函数(7)代入(5)，再代回到(4)式，利用式(8)做测试函数，最终的完全离散化时域方程为

3.2 无条件稳定(UCS)B-spline函数TDFEM

另一种选择是采用的B-spline函数和测试函数为

最终的完全离散化时域方程形如式(9)，参数 βi不同:β－1= β2=1/ 8，β0= β1=3/8.

3.3 稳定性分析

包含PML部分的TDFEM完全离散方程的稳定性分析很复杂，文献［13］中证明了式(4)中的［M］，{p}和{q}这几项，对整个方程的稳定性影响很微小，故而为了简便起见可以忽略不计.对化简了的式(4)，用Newmark-β法来离散化，再进行z变换，得到

定义 σ =max{σx，σy，σz}，式(10)可转化成特征值问题.假定矩阵［L］－1［S］Δt2的特征值为λ，式(10)可写成

稳定条件是式(11)中的根z落在单位圆内，相应的λ的取值范围在0～λmax，即

当β≥1/4时，不论Δt取多大都可以满足z的值在单位圆内，即该格式是无条件稳定的.

对式(9)引入β=2β2，用上面同样的方法对其进行z变换，得到的特征值方程与式(11)一样.所以，当在半离散方程(4)中采用Wn1(t)，时域离散方程(9)的参数β=13/60<1/ 4，是条件稳定的;而在半离散方程(4)中采用Wn2(t)，时域离散方程(9)的参数β=1/ 4，则是无条件稳定的.条件稳定方案的稳定条件是选取满足下式的Δt，ρ(·)代表谱半径.

4 算例分析

基于上述TDFEM编制FORTRAN程序并在P4计算机上实现，用此程序计算了一些算例来验证算法的精确性和有效性.脉冲函数采用高斯脉冲，由下式给出:

其中t0=6T，f0=1/T=2 946 MHz.图1是高斯脉冲的波形和波谱图.

图1 高斯脉冲的波形和波谱图

4.1 谐振腔

为了证实新算法的精确性和有效性，第一个算例是尺寸为0.072 m×0.050 m×0.072 m的无损谐振腔.采用高斯脉冲源，先计算了电场的时域响应，然后通过时间响应的快速傅立叶变换得到谐振腔的谐振频率.表1列出了各种时域基函数的TDFEM谐振腔的TE模谐振频率和对应的解析解.

从该表中可以看出，条件稳定和无条件稳定的B-spline函数TDFEM都和解析解吻合得很好，从而证明了前面提出的方法的精确性和有效性.在精度上，条件稳定二次B-spline算法显示出更小的相对误差.

表1 谐振腔的TE模谐振频率

4.2 偶极子天线

第二个算例是偶极子天线，由图2所示.天线的两极长度L皆为25 mm，馈源沿y轴.采用均匀剖分，未知量总数为10 560.图2给出了偶极子天线的瞬态响应，观察点在r=(－0.036x，0.036y，0.018z)m处.

图2 电场值的y方向分量

为了进一步证实远场结果的正确性，通过傅立叶变换获得了2.9 GHz下的E－面方向图，见图 3，并与文献［15］和传统的频域有限元算法的结果进行比较.其中，实线代表采用条件稳定B-spline时域基函数的TDFEM，圆圈代表采用传统TDFEM的结果，星星代表采用无条件稳定B－spline时域基函数，由此可以看出新时域有限元算法的数值结果与之很好地吻合.

图3 E－面方向图

在运算效率方面，传统的TDFEM，时间步长取dt=5 ps，运算时间为 649.6 s;若取 dt=10 ps，在420步发散;B－spline函数条件稳定方案，时间步长取dt=5 ps，运算时间为605.969 s;若取dt=10 ps，在400步发散;而无条件稳定的方案，时间步长取dt=10 ps，运算时间为455.437 s.从而证实了无条件稳定算法可以取较大的时间步长，节约了运算时间.

4.3 超宽带分形蝴蝶天线

第3个算例是超宽带分形蝴蝶天线，其结构如图4所示(左图是传统的蝴蝶天线，中间是S=1情况下的分形蝴蝶天线，右图是S=2情况下的分形蝴蝶天线).L=h=0.032 m，馈源在天线中心处，沿y轴方向.3种天线未知量总数分别为10 458，36 211和20 992.采用条件稳定B－spline函数，时间步长取dt=5 ps.

图4 分形蝴蝶天线的结构图

在图5中给出了这3种天线在频率增长情况下的E－面方向图，很容易发现在1.2 GHz时，3种天线的差别不大;而在3 GHz或是更高的6 GHz或8 GHz，传统的蝴蝶天线衰减更多.

通过以上算例，表明TDFEM可以直接在时域模拟超宽带天线的特性，而避免了传统频域方法在扫频上所消耗的计算时间.在超宽带天线和雷达设计方面可得到广泛的应用.

5 结论

文中推导了2种基于二次B-spline函数的时域有限元法计算格式，一种是条件稳定的方案，另一种是无条件稳定的方案.对这两种计算格式进行稳定性理论分析，得出该算法的稳定性条件.通过算例分析，得出条件稳定的方案在时间步长取值较大时会发散;而无条件稳定的方案在时间步长选取上有更大的范围，所以运算时间少于前者.数值结果证实了基于二次B-spline函数的时域有限元法的有效性.

图5 传统和分形蝴蝶天线在不同频率下的E－面方向图

［1］CANGELLARIS A C，LIN C C，MEI K.Point-matched time domain finite element methods for electromagnetic radiation and scattering［J］.IEEE Transactions on Antennas and Propagation， 1987，35(10):1160－1173.

［2］CANGELLARIS A C.Time-domain finite methods for electromagnetic wave propagation and scattering［J］.IEEE Transactions on Magnetics， 1991，27(5):3780－3785.

［3］LEE J F，LEE R，CANGELLARIS A C.Time-domain finite-element methods［J］.IEEE Transactions on Antennas and Propagation， 1997，45(3):430－442.

［4］GEDNEY S D，NAVSARIWALA U.An unconditionally stable finite element time-domain solution of the vector wave equation［J］.IEEE Microwave and Guided Wave Letters， 1995，5(10):332－334.

［5］BERENGER J P.A perfectly matched layer for the absorption of electromagnetic waves［J］.Journal of Computational Physics， 1994，114:185 －200.

［6］SCHOENBERG I J.Contribution to the problem of approximation of equidistant data by analytic function［J］.Quarterly of Applied Mathematics， 1946，4(1):45 － 99，4(2):112－141.

［7］WANG P，XIA M Y，JIN J M，et al.Time domain integral equation solvers using quadratic B-spline temporal basis functions［J］.Microwave and Optical Technology Letters， 2007，49(5):1154－1159.

［8］HOMSUP N.An B-spline based high-order finite difference timedomain schemes for the Maxwell equations［C］//Proceedings of Asia-Pacific Microwave Conference.Sydney:［s.n.］，2000:962－965.

［9］SACKS Z S，KINGSLAND D M，LEE R，et al.A perfectly matched anisotropic absorber for use as an absorbing boundary condition［J］.IEEE Transactions on Antennas and Propagation， 1995，43(12):1460－1463.

［10］LEE J F，SACKS Z.Whitney elements time domain(WETD)methods［J］.IEEE Transactions on Magnetics， 1995，31(3):1325－1329.

［11］GREVILLE T N E.Theory and Application of Spline Functions［M］.New York:Academic，1969.

［12］NEWMARK N.A method of computation for structural dynamics［J］.Journal of the Engineering Mechanics Division，. 1959，85:67 －94.

［13］JIAO D，JIN J M，MICHIELSSEN E，RILEY D.Time－domain finite-element simulation of three－dimensional scattering and radiation problems using perfectly matched layers［J］.IEEE Transactions on Antennas and Propagation， 2003，51(2):296－305.

［14］WU X，ZHOU L Z.Application of B-spline temporal basis function in time-domain finite element method for three-dimensional EM radiation problem［C］//Proceedings of 2008 Asia-Pacific Symposium on EMC＆19th InternationalZurich Symposium on Electromagnetic Compatibility.Singapore:［s.n.］，2008:774－777.

［15］LOU Z，JIN J M.Modeling and simulation of broadband antennas using the time-domain finite element method［J］.IEEE Transactions on Antennas and Propagation， 2005，53(12):4099 －4110.

Application of quadratic B-spline temporal basis function in time-domain finite element method

WU Xia1，2，ZHOU Le-zhu1

(1.School of Electronic Engineering ＆ Computer Science，Peking University，Beijing 100871，China，wuxiaeecs00@gmail.com;2.School of Information Science and Engineering，Fudan University，Shanghai 200433，China)

In order to be applied to electromagnetic(EM)radiation problems，a novel time-domain finite element method(TDFEM)on the basis of quadratic B-spline temporal basis functions is proposed.The formulation time-domain finite element methods is reviewed first.Then two TDFEM schemes based on quadratic B－spine temporal basis functions are presented:conditional stability and unconditional stability.They are validated via canonical cases.The efficiencies of the two schemes in precision level and run time are compared.The stability analysis of the proposed formulations is summarized and then verified through numerical results.

time-domain finite element method;quadratic B-spline temporal basis functions;electromagnetic radiation problems;unconditional stable;Ultra-wideband

O441

A

0367－6234(2010)05－0836－05

2008－12－22.

吴 霞 (1982—)，女，博士研究生;

周乐柱 (1944—)，男，教授，博士生导师.

(编辑 张 宏)