二次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—)，男，教授，博士生導師.

(編輯 張 宏)