吳 霞，周樂柱

(1.北京大學(xué)信息科學(xué)技術(shù)學(xué)院，北京 100871，wuxiaeecs00@gmail.com;2.復(fù)旦大學(xué)信息科學(xué)與工程學(xué)院，上海 200433)

二次B-Spline時域基函數(shù)的TDFEM的應(yīng)用

吳 霞1，2，周樂柱1

(1.北京大學(xué)信息科學(xué)技術(shù)學(xué)院，北京 100871，wuxiaeecs00@gmail.com;2.復(fù)旦大學(xué)信息科學(xué)與工程學(xué)院，上海 200433)

提出了時域有限元法(TDFEM)的一種新的基函數(shù)—二次B-spline時域基函數(shù).首先簡述了時域有限元法的原理和基本公式;然后提出了新型的基于B-spline函數(shù)的條件穩(wěn)定和無條件穩(wěn)定的時域有限元法方案，并應(yīng)用于三維電磁輻射問題.通過典型的算例對這兩種方案的精度、運算時間進行了比較，證實了基于二次B-spline函數(shù)的時域有限元法的有效性.通過穩(wěn)定性理論分析得出該算法的精確穩(wěn)定性，并且通過數(shù)值計算的結(jié)果得到驗證.

時域有限元法;二次B-spline時域基函數(shù);電磁輻射;無條件穩(wěn)定;超寬帶

在計算電磁領(lǐng)域，有不同的時域分析計算方法，例如時域有限差分法(FDTD)、時域積分法(TDIE)和時域有限元法(TDFEM)等，致力于對各種不同的結(jié)構(gòu)進行建模與仿真分析.在眾多的方法中，TDFEM因其繼承了FEM特別適用于復(fù)雜的幾何結(jié)構(gòu)和介電特性分布的優(yōu)點，及在非均勻介質(zhì)體內(nèi)產(chǎn)生的也是稀疏矩陣的長處，因而在解決寬帶特性的電磁場問題上彌補了FDTD和TDIE法的某些不足而被人們漸漸關(guān)注.TDFEM起源于上世紀(jì)80年代，分為兩大類，一類是最初由 A.C.Cangellaris等［1－3］提出的微分形式的點匹配時域有限元法.該方法結(jié)合了FDTD顯式積分法的簡明性和傳統(tǒng)有限元的靈活性，但是會出現(xiàn)源于電場和磁場點在網(wǎng)格上的交替出現(xiàn)而產(chǎn)生的類似FDTD法的越級(leap-frog)問題.為了避免或解決這些問題，第2類隱式時域有限元法在1995年提出［4］.該方法是基于電場(或磁場)的矢量有限元空間基函數(shù)和時間基函數(shù)展開的二階矢量波方程，其優(yōu)點是簡單而穩(wěn)定，缺點是每個時間步必須求解一個矩陣方程;對于開域空間問題，采用傳統(tǒng)的吸收邊界條件，有限元計算空間十分龐大，這個問題更加嚴(yán)重，以致無法進行.由于這些原因，這類時域有限元法發(fā)展緩慢.直到2000年初，由于完全匹配層(PML)［5］等概念的提出和應(yīng)用，時域有限元法才得到較大的發(fā)展.

本文研究上述第2類的時域有限元法，與已報道的TDFEM方法不同的是，本文首次引入了B-spline函數(shù)作為時間基函數(shù).另外，利用Newmark-β方法驗證了相應(yīng)的檢測函數(shù).B-spline函數(shù)早在1946年被I.J.Schenberg提出［6］，因其插值特性和計算簡單的優(yōu)點近期才在計算電磁算法中被重視.B-spline時域基函數(shù)已經(jīng)在FDTD和TDIE算法中被采用并取得了不錯的效果［7－8］.然而，目前在時域有限元法的研究中還沒有具體應(yīng)用B-spline時域基函數(shù)的報道.本文算法的另外兩個特點是:1)采用Sacks PML做截斷邊界條件;2)采用高階矢量有限元作為空間基函數(shù).

1 TDFEM基本公式

為了推導(dǎo)出電場的時域有限元方法，從一個時域電磁場的初值問題開始［3］，方程式如下:

其中J0代表外加電流源，εr和μr分別是相對介電常數(shù)和相對磁導(dǎo)率，c0代表自由空間光速.用Sacks PML來截斷有限區(qū)域，在PML區(qū)域，電場滿足下式:

為了把式(4)的半離散化的時域方程轉(zhuǎn)變成完全離散化的時域方程，ej(t)應(yīng)當(dāng)由某種時間基函數(shù)Tm(t)展開，即

其中emj是ej(t)在t=tm(m為整數(shù))時相對應(yīng)的值.在以往報道的時域有限元法文獻中，Tm(t)都采用如下的線性時域基函數(shù):

本文采用二次B－spline函數(shù)作為時域基函數(shù).

2 B-spline時域基函數(shù)

構(gòu)造B樣條函數(shù)時，給定階數(shù)n，B樣條函數(shù)φn(x)的表達式為［11］

其中:H(x)是階躍函數(shù).B樣條函數(shù)雖然是分段函數(shù)，但(n－2)次導(dǎo)數(shù)連續(xù)，因此用在TDFEM中可以更精確地進行時域離散.

在式(6)中取n=3得到本文所用的二階B-spline時域基函數(shù)的表達式由下式給出.

3 兩種基于B-spline的TDFEM

該新算法的穩(wěn)定性是至關(guān)重要的.鑒于穩(wěn)定性依賴于時域的離散化，而且有多種方法把式(4)變成一個完全離散的矩陣方程，如中心差分、前向差分、后向差分或Newmark-β 法［12］.本文采取和空間離散同樣的方法，即用時域基函數(shù)來獲得完全離散化的方程.時域離散展開系數(shù)的值取決于離散化時域基函數(shù)和測試函數(shù)的具體形式.

3.1 條件穩(wěn)定(CS)B-spline函數(shù)TDFEM方案

第一種選擇是采用如下的B-spline函數(shù)和測試函數(shù):

將時域基函數(shù)(7)代入(5)，再代回到(4)式，利用式(8)做測試函數(shù)，最終的完全離散化時域方程為

3.2 無條件穩(wěn)定(UCS)B-spline函數(shù)TDFEM

另一種選擇是采用的B-spline函數(shù)和測試函數(shù)為

最終的完全離散化時域方程形如式(9)，參數(shù) βi不同:β－1= β2=1/ 8，β0= β1=3/8.

3.3 穩(wěn)定性分析

包含PML部分的TDFEM完全離散方程的穩(wěn)定性分析很復(fù)雜，文獻［13］中證明了式(4)中的［M］，{p}和{q}這幾項，對整個方程的穩(wěn)定性影響很微小，故而為了簡便起見可以忽略不計.對化簡了的式(4)，用Newmark-β法來離散化，再進行z變換，得到

定義 σ =max{σx，σy，σz}，式(10)可轉(zhuǎn)化成特征值問題.假定矩陣［L］－1［S］Δt2的特征值為λ，式(10)可寫成

穩(wěn)定條件是式(11)中的根z落在單位圓內(nèi)，相應(yīng)的λ的取值范圍在0～λmax，即

當(dāng)β≥1/4時，不論Δt取多大都可以滿足z的值在單位圓內(nèi)，即該格式是無條件穩(wěn)定的.

對式(9)引入β=2β2，用上面同樣的方法對其進行z變換，得到的特征值方程與式(11)一樣.所以，當(dāng)在半離散方程(4)中采用Wn1(t)，時域離散方程(9)的參數(shù)β=13/60<1/ 4，是條件穩(wěn)定的;而在半離散方程(4)中采用Wn2(t)，時域離散方程(9)的參數(shù)β=1/ 4，則是無條件穩(wěn)定的.條件穩(wěn)定方案的穩(wěn)定條件是選取滿足下式的Δt，ρ(·)代表譜半徑.

4 算例分析

基于上述TDFEM編制FORTRAN程序并在P4計算機上實現(xiàn)，用此程序計算了一些算例來驗證算法的精確性和有效性.脈沖函數(shù)采用高斯脈沖，由下式給出:

其中t0=6T，f0=1/T=2 946 MHz.圖1是高斯脈沖的波形和波譜圖.

圖1 高斯脈沖的波形和波譜圖

4.1 諧振腔

為了證實新算法的精確性和有效性，第一個算例是尺寸為0.072 m×0.050 m×0.072 m的無損諧振腔.采用高斯脈沖源，先計算了電場的時域響應(yīng)，然后通過時間響應(yīng)的快速傅立葉變換得到諧振腔的諧振頻率.表1列出了各種時域基函數(shù)的TDFEM諧振腔的TE模諧振頻率和對應(yīng)的解析解.

從該表中可以看出，條件穩(wěn)定和無條件穩(wěn)定的B-spline函數(shù)TDFEM都和解析解吻合得很好，從而證明了前面提出的方法的精確性和有效性.在精度上，條件穩(wěn)定二次B-spline算法顯示出更小的相對誤差.

表1 諧振腔的TE模諧振頻率

4.2 偶極子天線

第二個算例是偶極子天線，由圖2所示.天線的兩極長度L皆為25 mm，饋源沿y軸.采用均勻剖分，未知量總數(shù)為10 560.圖2給出了偶極子天線的瞬態(tài)響應(yīng)，觀察點在r=(－0.036x，0.036y，0.018z)m處.

圖2 電場值的y方向分量

為了進一步證實遠場結(jié)果的正確性，通過傅立葉變換獲得了2.9 GHz下的E－面方向圖，見圖 3，并與文獻［15］和傳統(tǒng)的頻域有限元算法的結(jié)果進行比較.其中，實線代表采用條件穩(wěn)定B-spline時域基函數(shù)的TDFEM，圓圈代表采用傳統(tǒng)TDFEM的結(jié)果，星星代表采用無條件穩(wěn)定B－spline時域基函數(shù)，由此可以看出新時域有限元算法的數(shù)值結(jié)果與之很好地吻合.

圖3 E－面方向圖

在運算效率方面，傳統(tǒng)的TDFEM，時間步長取dt=5 ps，運算時間為 649.6 s;若取 dt=10 ps，在420步發(fā)散;B－spline函數(shù)條件穩(wěn)定方案，時間步長取dt=5 ps，運算時間為605.969 s;若取dt=10 ps，在400步發(fā)散;而無條件穩(wěn)定的方案，時間步長取dt=10 ps，運算時間為455.437 s.從而證實了無條件穩(wěn)定算法可以取較大的時間步長，節(jié)約了運算時間.

4.3 超寬帶分形蝴蝶天線

第3個算例是超寬帶分形蝴蝶天線，其結(jié)構(gòu)如圖4所示(左圖是傳統(tǒng)的蝴蝶天線，中間是S=1情況下的分形蝴蝶天線，右圖是S=2情況下的分形蝴蝶天線).L=h=0.032 m，饋源在天線中心處，沿y軸方向.3種天線未知量總數(shù)分別為10 458，36 211和20 992.采用條件穩(wěn)定B－spline函數(shù)，時間步長取dt=5 ps.

圖4 分形蝴蝶天線的結(jié)構(gòu)圖

在圖5中給出了這3種天線在頻率增長情況下的E－面方向圖，很容易發(fā)現(xiàn)在1.2 GHz時，3種天線的差別不大;而在3 GHz或是更高的6 GHz或8 GHz，傳統(tǒng)的蝴蝶天線衰減更多.

通過以上算例，表明TDFEM可以直接在時域模擬超寬帶天線的特性，而避免了傳統(tǒng)頻域方法在掃頻上所消耗的計算時間.在超寬帶天線和雷達設(shè)計方面可得到廣泛的應(yīng)用.

5 結(jié)論

文中推導(dǎo)了2種基于二次B-spline函數(shù)的時域有限元法計算格式，一種是條件穩(wěn)定的方案，另一種是無條件穩(wěn)定的方案.對這兩種計算格式進行穩(wěn)定性理論分析，得出該算法的穩(wěn)定性條件.通過算例分析，得出條件穩(wěn)定的方案在時間步長取值較大時會發(fā)散;而無條件穩(wěn)定的方案在時間步長選取上有更大的范圍，所以運算時間少于前者.數(shù)值結(jié)果證實了基于二次B-spline函數(shù)的時域有限元法的有效性.

圖5 傳統(tǒng)和分形蝴蝶天線在不同頻率下的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，F(xiàn)udan 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—)，男，教授，博士生導(dǎo)師.

(編輯 張 宏)