吳 霞，周樂柱

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

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

吳 霞1，2，周樂柱1

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

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

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

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

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

1 TDFEM基本公式

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

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

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

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

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

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

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

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

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

3 兩種基于B-spline的TDFEM

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

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

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

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

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

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

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

3.3 穩(wěn)定性分析

包含PML部分的TDFEM完全離散方程的穩(wěn)定性分析很復(fù)雜，文獻(xiàn)［13］中證明了式(4)中的［M］，{p}和{q}這幾項(xiàng)，對(duì)整個(gè)方程的穩(wěn)定性影響很微小，故而為了簡便起見可以忽略不計(jì).對(duì)化簡了的式(4)，用Newmark-β法來離散化，再進(jìn)行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時(shí)，不論Δt取多大都可以滿足z的值在單位圓內(nèi)，即該格式是無條件穩(wěn)定的.

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

4 算例分析

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

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

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

4.1 諧振腔

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

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

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

4.2 偶極子天線

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

圖2 電場值的y方向分量

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

圖3 E－面方向圖

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

4.3 超寬帶分形蝴蝶天線

第3個(gè)算例是超寬帶分形蝴蝶天線，其結(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ù)，時(shí)間步長取dt=5 ps.

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

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

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

5 結(jié)論

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

(編輯 張 宏)