利用時域間斷伽遼金算法求解Tellegen介質中的電磁波

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

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

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

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

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

1 微分方程公式

1.1 麥克斯韋方程組

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

(1)

(2)

Tellegen介質的本構關系方程為：

D=εE+χH

(3)

B=χE+μH

(4)

式中：χ為Tellegen參數。該參數滿足條件：

χ2<με

(5)

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

(7)

(8)

1.2 Tellegen介質的弱解形式

根據以上公式,我們推導出適用于Tellegen介質的麥克斯韋方程的弱解形式,如下：

(9)

(10)

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

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

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

(13)

式中：M為自由度的個數。

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

1.3 Tellegen介質的系統矩陣方程

進一步將半離散化方程組總結為系統矩陣方程組：

其中,系統矩陣的表達式為：

(18)

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

1.4 時間步進方案

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

(20)

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

(22)

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

1.5 偏轉角度的計算

根據文獻[20],我們可以推導出在Tellegen介質的半空間模型中的反射系數與透射系數[20],其表達式為：

因此,反射偏轉角度與透射偏轉角度的計算式為：

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

2 數值算例

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

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

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

圖2 計算域的立方體網格剖分圖

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

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

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

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

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

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

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

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

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

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

表1 DGTD模型觀測點處的偏轉角度與文獻[1]中的結果對比

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

(a) 側視圖

(a) 側視圖

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

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

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

表3 DGTD模型觀測點處的偏轉角度與解析解對比

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

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

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

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

3 結語

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