利用時(shí)域間斷伽遼金算法求解Tellegen介質(zhì)中的電磁波

黃 鵬,周遠(yuǎn)國(guó),* 任 強(qiáng),王 焱

(1.西安科技大學(xué)通信與信息工程學(xué)院,西安,710054;2.北京航空航天大學(xué)電子信息工程學(xué)院,北京,100191;3.中國(guó)航空研究院研究生院,江蘇揚(yáng)州,226000;4.電磁環(huán)境效應(yīng)航空科技重點(diǎn)實(shí)驗(yàn)室,沈陽(yáng),110035)

在宏觀(guān)層面,雙各向同性介質(zhì)的本構(gòu)關(guān)系由兩個(gè)磁電耦合的參數(shù)表征。具體來(lái)說(shuō),雙各向同性介質(zhì)又分為手性介質(zhì)和Tellegen介質(zhì)[1]。手性介質(zhì)具有互異性和固有的手性特征,以及眾所周知的圓雙折射性和二色性現(xiàn)象。而Tellegen介質(zhì)更加復(fù)雜[2][3],呈現(xiàn)非互易性,它與激發(fā)場(chǎng)呈同相磁電耦合特征。在普通介質(zhì)中,電通量密度D僅與電場(chǎng)強(qiáng)度E有關(guān),磁通量密度B僅與磁場(chǎng)強(qiáng)度H有關(guān);然而,在Tellegen介質(zhì)中,它的電通量密度D不僅與電場(chǎng)強(qiáng)度E有關(guān),還與磁場(chǎng)強(qiáng)度H有關(guān);同樣地,磁通量密度也是如此。

另一方面,時(shí)域間斷伽遼金算法(time-domain discontinuous Galerkin method,DGTD)作為一種高效的數(shù)值方法,在計(jì)算電磁學(xué)中快速興起[4-5]。該技術(shù)利用伽遼金加權(quán)法離散麥克斯韋方程組的弱解形式,放松了單元之間的邊界條件,單元之間通過(guò)數(shù)值通量相互聯(lián)系并交換數(shù)據(jù)[6]?？刹捎萌切巍⑺倪呅巍⑺拿骟w或六面體作為網(wǎng)格單元[7],并可以靈活選擇結(jié)點(diǎn)基函數(shù)或者棱邊基函數(shù)構(gòu)建系統(tǒng)矩陣。需要高精度時(shí)采用高階基函數(shù),對(duì)精度要求低時(shí)采用低階基函數(shù)。這些特點(diǎn)使得DGTD可以靈活的分析復(fù)雜幾何問(wèn)題,復(fù)雜介質(zhì)問(wèn)題以及多尺度問(wèn)題。據(jù)研究,2013年格拉納達(dá)大學(xué)的J.A.Gonzalez博士系統(tǒng)闡述了DGTD理論的發(fā)展,從Maxwell方程組到DGTD時(shí)域步進(jìn)公式[8],并詳細(xì)討論了黎曼邊界條件導(dǎo)出的數(shù)值通量,分析了不同階基函數(shù)情況下的算法精確度[9]。2015至2016年,Ren等[12]采用了區(qū)域分解技術(shù)解決了電大尺寸目標(biāo)計(jì)算資源消耗過(guò)大的問(wèn)題[10-11]。同年,楊謙等人利用時(shí)域間斷伽遼金算法計(jì)算了空腔與填充諧振腔諧振頻率[12]。2017年Su[13]等提出了基于動(dòng)態(tài)自適應(yīng)笛卡爾網(wǎng)格的DGTD方法,用于解決電磁場(chǎng)的全波仿真問(wèn)題[13]。2018年,買(mǎi)文鼎等人提出了一種改進(jìn)的自適應(yīng)判據(jù),可以兼顧精度與效率,有效地處理電磁全波計(jì)算問(wèn)題[14]。2021年Li等[15]人評(píng)述并討論了多物理 DGTD仿真問(wèn)題,驗(yàn)證了 DGTD方法對(duì)多物理場(chǎng)耦合問(wèn)題的用適性[16]。2022年,Chen L等[17]首次推導(dǎo)了廣義表面?zhèn)鬏敆l件模型(generalized sheet transition conditions,GSTCs)的數(shù)值通量,并用于DGTD算法,從而有效地提高了超表面問(wèn)題的仿真效率[18]。

然而,迄今為止國(guó)內(nèi)外對(duì)Tellegen介質(zhì)中電磁波傳播問(wèn)題研究甚少。僅文獻(xiàn)[1]記錄2014年西班牙A.Grande博士利用時(shí)域有限差分方法(finite difference time domain,FDTD)對(duì)Tellegen介質(zhì)進(jìn)行了研究,分析了平面波在這種介質(zhì)中的極化偏轉(zhuǎn)角度。FDTD算法是一種發(fā)展較為成熟,且應(yīng)用廣泛的數(shù)值算法,例如利用該方法分析共面波導(dǎo)傳輸結(jié)構(gòu)[19]、以及隱身目標(biāo)超寬帶雙站RCS的計(jì)算等[20]。然而該算法也有其固有的缺點(diǎn),即它的階梯型誤差和數(shù)值色散的問(wèn)題。相比于該算法,DGTD具有網(wǎng)格靈活、高階精度、適合多尺度計(jì)算等優(yōu)點(diǎn),并且有效克服了FDTD方法的缺點(diǎn)。為了處理雙各向同性介質(zhì)中大規(guī)模的多尺度問(wèn)題,本文針對(duì)Tellegen介質(zhì)的本構(gòu)關(guān)系,推導(dǎo)了相應(yīng)的麥克斯韋方程弱解形式,提出了一種新型的DGTD系統(tǒng)矩陣的離散方案,準(zhǔn)確模擬了平面波在Tellegen介質(zhì)中的瞬態(tài)傳播特性。這項(xiàng)工作為仿真雙各向同性介質(zhì)中大規(guī)模地多尺度電磁問(wèn)題打下基礎(chǔ)。

1 微分方程公式

1.1 麥克斯韋方程組

首先,從麥克斯韋方程組的微分形式出發(fā)：

(1)

(2)

Tellegen介質(zhì)的本構(gòu)關(guān)系方程為：

D=εE+χH

(3)

B=χE+μH

(4)

式中：χ為T(mén)ellegen參數(shù)。該參數(shù)滿(mǎn)足條件：

χ2<με

(5)

式中：χr為相對(duì)Tellegen參數(shù)。將式(3)～(4)代入式(1)、式(2),可以得到：

(7)

(8)

1.2 Tellegen介質(zhì)的弱解形式

根據(jù)以上公式,我們推導(dǎo)出適用于Tellegen介質(zhì)的麥克斯韋方程的弱解形式,如下：

(9)

(10)

利用分部積分法,進(jìn)一步得到：

式中：νe、νh分別表示電場(chǎng)試探函數(shù)、磁場(chǎng)試探函數(shù)。

假設(shè)電場(chǎng)E與磁場(chǎng)H可由矢量基函數(shù)Φ表示為如下形式：

(13)

式中：M為自由度的個(gè)數(shù)。

結(jié)合式(11)～式(13),通過(guò)間斷伽遼金半離散化可獲得如下方程：

1.3 Tellegen介質(zhì)的系統(tǒng)矩陣方程

進(jìn)一步將半離散化方程組總結(jié)為系統(tǒng)矩陣方程組：

其中,系統(tǒng)矩陣的表達(dá)式為：

(18)

式中：Zi、Yi分別表示波阻抗與導(dǎo)納;ni為單元表面上指向朝外的法向量。在上述的離散化系統(tǒng)矩陣方程中,單元與單元之間進(jìn)行信息交換的數(shù)值通量采用的是迎風(fēng)通量[19]。

1.4 時(shí)間步進(jìn)方案

本文采用龍格庫(kù)塔法,對(duì)式(16)～式(17)所示的系統(tǒng)矩陣方程組進(jìn)行迭代求解。四階龍格庫(kù)塔公式如下所示：

(20)

由上述公式,基于四階龍格庫(kù)塔公式的時(shí)間步進(jìn)方案如下所示：

(22)

式中：ak,l、bk、ck均為 Butcher 表中的系數(shù)。式(21)中的矩陣表達(dá)式為：

1.5 偏轉(zhuǎn)角度的計(jì)算

根據(jù)文獻(xiàn)[20],我們可以推導(dǎo)出在Tellegen介質(zhì)的半空間模型中的反射系數(shù)與透射系數(shù)[20],其表達(dá)式為：

因此,反射偏轉(zhuǎn)角度與透射偏轉(zhuǎn)角度的計(jì)算式為：

此外,將θEH定義為平面波的電場(chǎng)極化方向與磁場(chǎng)方向的夾角,它滿(mǎn)足如下關(guān)系式：

2 數(shù)值算例

為了模擬Tellegen介質(zhì)中平面波的傳播特性,首先考慮Air-Tellegen分層空間模型,如圖1所示。平面波由空氣垂直入射到Tellegen介質(zhì)中,并假設(shè)該分層介質(zhì)的邊界為無(wú)窮遠(yuǎn)。

圖1 包含空氣與Tellegen介質(zhì)的分層空間模型

該模型的網(wǎng)格剖分圖如圖2所示,計(jì)算域的大小為0.018 5 m×0.05 m×0.05 m,其中綠色部分為空氣介質(zhì)的網(wǎng)格剖分,黃色部分為T(mén)ellegen介質(zhì)的網(wǎng)格剖分,2種介質(zhì)區(qū)域內(nèi)均采用立方體網(wǎng)格進(jìn)行剖分。值得注意的是,2種區(qū)域內(nèi)系統(tǒng)矩陣的運(yùn)算是獨(dú)立的,空氣與Tellegen介質(zhì)之間的交界面處僅通過(guò)迎風(fēng)通量來(lái)實(shí)現(xiàn)區(qū)域之間的信息交互。平面波采用余弦調(diào)制高斯脈沖信號(hào)作為激勵(lì)源沿x軸方向垂直入射,電場(chǎng)極化方向?yàn)閦軸正方向,中心頻率為 9 GHz。Tellegen介質(zhì)的電磁參數(shù)分別為εr=3.5,μr=1.2,χr=0.6。在計(jì)算區(qū)域內(nèi)設(shè)置2個(gè)觀(guān)測(cè)點(diǎn)：觀(guān)測(cè)點(diǎn)1位于空氣中,距離分界面 0.09 m;觀(guān)測(cè)點(diǎn)2位于Tellegen介質(zhì)中,距離分界面 -0.012 5 m。

圖2 計(jì)算域的立方體網(wǎng)格剖分圖

利用DGTD進(jìn)行仿真時(shí),電場(chǎng)E的基函數(shù)階數(shù)為三階,磁場(chǎng)H的基函數(shù)階數(shù)為二階,時(shí)間步長(zhǎng)dt為0.000 5×10-9s,時(shí)間迭代的步數(shù)為4 000。數(shù)值實(shí)驗(yàn)分別記錄了觀(guān)測(cè)點(diǎn)1和2的電場(chǎng)分量變化情況,如圖3、圖4所示。

圖3 觀(guān)測(cè)點(diǎn)1處的電場(chǎng)變化

圖4 觀(guān)測(cè)點(diǎn)2處的電場(chǎng)變化

從圖3可以看到,由于z方向極化的原因,觀(guān)測(cè)點(diǎn)1的直達(dá)波只有電場(chǎng)z分量,而在反射波中出現(xiàn)了y分量。也就是說(shuō),當(dāng)平面波垂直入射到Tellegen介質(zhì)時(shí),它會(huì)改變反射波的電場(chǎng)極化方向使其發(fā)生偏轉(zhuǎn)。同樣地,Tellegen介質(zhì)也會(huì)對(duì)透射波的極化方向產(chǎn)生影響,圖4可以看出透射波中包含電場(chǎng)y方向分量。

為了進(jìn)一步觀(guān)察反射波和透射波的偏轉(zhuǎn)情況,分別繪制了觀(guān)測(cè)點(diǎn)1和觀(guān)測(cè)點(diǎn)2的三維電場(chǎng)分量圖,如圖5和6所示。

圖6 觀(guān)測(cè)點(diǎn)2處的電場(chǎng)變化值(三維)

從圖中可以看到,由于Tellegen介質(zhì)本構(gòu)關(guān)系中電磁場(chǎng)交叉耦合,使得照射到該介質(zhì)上的電磁波偏振方向發(fā)生了明顯的改變。

圖7 (a)為觀(guān)測(cè)點(diǎn)1處入射波與反射波的極化偏轉(zhuǎn)情況。圖7 (b)為觀(guān)測(cè)點(diǎn)2處入射波與透射波的極化偏轉(zhuǎn)情況。通過(guò)這兩幅圖,可以更直觀(guān)地觀(guān)察Tellegen介質(zhì)對(duì)平面波極化方向的調(diào)制。

(a)觀(guān)測(cè)點(diǎn)1(正視圖) b) 觀(guān)測(cè)點(diǎn)2(正視圖)

為了定量研究Tellegen介質(zhì)對(duì)電磁波極化方向的影響,利用DGTD算法計(jì)算了反射波與透射波的極化偏轉(zhuǎn)角度,并與文獻(xiàn)[1]中的結(jié)果進(jìn)行了對(duì)比。從表1可以看到DGTD計(jì)算反射波與透射波的極化偏轉(zhuǎn)角分別為27.771°和-10.388°,與參考文獻(xiàn)[1]中結(jié)果吻合較好。

表1 DGTD模型觀(guān)測(cè)點(diǎn)處的偏轉(zhuǎn)角度與文獻(xiàn)[1]中的結(jié)果對(duì)比

第2個(gè)模型算例,也為Air-Tellegen分層空間模型,唯一不同的是這個(gè)算例中的電磁參數(shù)為εr=3,μr=2.1,χr=0.5。其模型的大小、基函數(shù)階數(shù)、時(shí)間步長(zhǎng)dt和時(shí)間迭代步數(shù)等參數(shù)與算例1相同。平面波的中心頻率為9 GHz,模型中的觀(guān)測(cè)點(diǎn)分別位于(0.09,0,0)、(-0.002 5,0,0),然后利用DGTD程序進(jìn)行仿真,其仿真結(jié)果圖8、圖9所示。

(a) 側(cè)視圖

(a) 側(cè)視圖

從圖8、圖9可以看出,算例2中的Tellegen介質(zhì)對(duì)平面波的影響,具有類(lèi)似的效果,只是反射偏轉(zhuǎn)角度與透射偏轉(zhuǎn)角度大小不同,通過(guò)計(jì)算電場(chǎng)極化方向與磁場(chǎng)方向之間的夾角,即θEH,算例2的計(jì)算結(jié)果與理論值的比較如表2所示。

表2 DGTD模型觀(guān)測(cè)點(diǎn)處的偏轉(zhuǎn)角度與理論值之間的對(duì)比

電場(chǎng)極化偏轉(zhuǎn)角φeT,磁場(chǎng)偏轉(zhuǎn)角φhT,以及反射偏轉(zhuǎn)角φr,得到的結(jié)果與解析解對(duì)比,如表3所示,可以看到DGTD的計(jì)算結(jié)果與理論值吻合較好。

表3 DGTD模型觀(guān)測(cè)點(diǎn)處的偏轉(zhuǎn)角度與解析解對(duì)比

第3個(gè)算例為 “彈頭”模型,彈頭的電磁參數(shù)為：εr=3.5,μr=1.2,χr=0.6,如圖10所示。入射波的方向從右向左,即(-1,0,0)。入射波的中心頻率為155 MHz,激勵(lì)脈沖采用余弦調(diào)制的高斯脈沖信號(hào)。分別在(2.5 1.3 1.3)、(-2.5 0 0)處設(shè)置觀(guān)測(cè)點(diǎn),觀(guān)察彈頭內(nèi)外的電場(chǎng)變化情況。為了更好地模擬彈頭的曲面結(jié)構(gòu),上述算例模型采用四面體網(wǎng)格剖分。

利用DGTD計(jì)算時(shí),時(shí)間步長(zhǎng)設(shè)置為dt= 0.05×10-9,時(shí)間步數(shù)nt= 3 000。計(jì)算可得觀(guān)測(cè)點(diǎn)處的電場(chǎng)變化情況如圖11所示。

(a) 觀(guān)測(cè)點(diǎn)1(2.5,1.3,1.3)處的電場(chǎng)變化情況

由圖11的結(jié)果可以看出,在彈頭模型外的觀(guān)測(cè)點(diǎn)處,即接收點(diǎn)1,首先記錄到入射場(chǎng),然后記錄由彈頭表面所形成的各級(jí)反射波,如圖11 (a)所示;在彈頭內(nèi)的觀(guān)測(cè)點(diǎn),即接收點(diǎn)2,首先記錄進(jìn)入Tellegen介質(zhì)中的透射波。由于電磁場(chǎng)在Tellegen區(qū)域內(nèi)的多重散射現(xiàn)象,接收點(diǎn)2記錄的波形振蕩比較嚴(yán)重,如圖11 (b)所示。

3 結(jié)語(yǔ)

基于時(shí)域間斷伽遼金(DGTD)離散方法,本文首次推導(dǎo)了適用Tellegen介質(zhì)的系統(tǒng)矩陣方程,并采用迎風(fēng)通量實(shí)現(xiàn)計(jì)算域中各個(gè)單元之間的信息交互。通過(guò)數(shù)值算例,我們計(jì)算并分析了平面波在具有雙電磁耦合效應(yīng)的Tellegen介質(zhì)中的傳播特性,驗(yàn)證了該算法的可行性和有效性。