黃 鵬,周遠國,* 任 強,王 焱

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

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

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

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

1 微分方程公式

1.1 麥克斯韋方程組

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

(1)

(2)

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

D=εE+χH

(3)

B=χE+μH

(4)

式中：χ為Tellegen參數(shù)。該參數(shù)滿足條件：

χ2<με

(5)

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

(7)

(8)

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

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

(9)

(10)

利用分部積分法,進一步得到：

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

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

(13)

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

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

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

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

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

(18)

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

1.4 時間步進方案

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

(20)

由上述公式,基于四階龍格庫塔公式的時間步進方案如下所示：

(22)

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

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

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

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

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

2 數(shù)值算例

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

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

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

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

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

圖3 觀測點1處的電場變化

圖4 觀測點2處的電場變化

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

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

圖6 觀測點2處的電場變化值(三維)

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

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

(a)觀測點1(正視圖) b) 觀測點2(正視圖)

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

表1 DGTD模型觀測點處的偏轉(zhuǎn)角度與文獻[1]中的結(jié)果對比

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

(a) 側(cè)視圖

(a) 側(cè)視圖

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

表2 DGTD模型觀測點處的偏轉(zhuǎn)角度與理論值之間的對比

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

表3 DGTD模型觀測點處的偏轉(zhuǎn)角度與解析解對比

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

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

(a) 觀測點1(2.5,1.3,1.3)處的電場變化情況

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

3 結(jié)語

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