適用于一維精細結構電磁目標模擬的通用HIE-FDTD 方法及程序實現*

牟春暉 陳娟2)? 范凱航 魯藝

1) (西安交通大學信息與通信工程學院,西安 710049)

2) (西安交通大學深圳研究院,深圳 518057)

混合顯-隱式時域有限差分方法的時間步長只由空間兩個方向上的網格大小決定,因而該方法被廣泛應用于沿一個方向具有精細結構的電磁目標的模擬.本文通過對傳統時域有限差分方法基本公式進行近似,給出適用于線性、非色散空間,一維精細結構電磁目標模擬的通用混合顯-隱式時域有限差分方法.在通用混合顯-隱式時域有限差分方法中,通過媒質編號索引和迭代系數空間定位實現三對角矩陣的自動構建,引入等效參數方法解決麥克斯韋方程微分形式在不同介質分界面失效的問題,引入卷積完全匹配層邊界條件對計算區域進行截斷.此外,本文提出了一種縮減三對角矩陣數量的方法,該方法可有效降低程序運行所需內存,提高程序運行效率.應用通用混合顯-隱式時域有限差分方法仿真了平面波照射介質板和雙頻微帶倒F 天線兩個模型,數值計算結果與傳統時域有限差分方法和CST 軟件的計算結果一致,但與傳統時域有限差分方法相比,通用混合顯-隱式時域有限差分方法的計算效率大大提高.

1 引言

時域有限差分(finite-difference time-domain,FDTD)方法[1,2]自提出以來已被廣泛應用于電磁學研究的眾多領域[3-9].傳統FDTD 方法因采用顯式差分算法而必須滿足Courant-Friedrich-Levy時間穩定性條件[2].因此,傳統FDTD 方法的時間步長受制于計算空間網格的最小值.當模擬目標在一個方向上具有精細結構時,為了保證計算精度,在目標精細結構所在的方向需要采用精細網格,依據時間穩定性條件,模擬此類問題時,時間步長需要取得很小,這就使得總的計算步數劇增.混合顯-隱式FDTD(hybrid implicit-explicit FDTD,HIEFDTD)方法[10-12]通過對精細網格所在方向的空間求導項采用混合顯-隱式差分算法,從而使得時間步長不再受限于精細網格大小.HIE-FDTD 方法自提出以來已被廣泛應用于模擬在一個方向具有精細結構的問題,如基于石墨烯的各種電磁器件[13-16]、微帶電路和天線[17,18]、線性/非線性集總網絡[19,20]等.

在使用FDTD 方法計算空間任意一個網格上的電場和磁場時,需要用到該網格處介質的相關參數.因此在FDTD 方法的基本公式中,所有的迭代系數均需要標記空間位置.此外,由于微分形式的麥克斯韋方程在介質參數突變面失效,因此在編寫FDTD 程序時需要引入等效參數[21,22]概念.然而,現有的關于HIE-FDTD 方法的文獻和專著中所給出的HIE-FDTD 基本公式,大多數沒有對迭代系數進行空間位置標記.依據這些公式編寫的HIE-FDTD 程序,在仿真不同的電磁目標時,需根據目標介質分布情況手動構建三對角矩陣,這就導致HIE-FDTD 程序不具備通用性.少數文獻在HIE-FDTD 迭代系數中增加了空間位置標記,但沒有考慮介質分界面處微分形式麥克斯韋方程失效問題,在計算過程中也沒有設置吸收邊界,因此在模擬復雜目標時,計算精度難以保證,且不能計算開放區域問題[19,20].

本文提出了一種適用于線性、非色散空間,一維精細結構電磁目標模擬的通用HIE-FDTD 程序實現方法.在HIE-FDTD 基本公式中,對所有迭代系數進行空間位置標記,并對模擬空間中的所有媒質進行編號.在介質分界面引入等效參數概念.同時,引入卷積完全匹配層(convolutional perfectly matched layer,CPML)邊界條件來截斷計算區域.依據通用HIE-FDTD 方法編寫的程序可根據目標介質分布情況自動構建三對角矩陣,且既可以模擬封閉區域問題,又可以模擬開放區域問題.此外,提出了一種縮減三對角矩陣數量的方法,該方法可有效降低程序運行所需內存,提高程序運行效率.本文應用通用HIE-FDTD 程序分別計算了平面波照射介質板時的透射場和反射場以及雙頻微帶倒F天線S參數和輻射方向圖,仿真結果與傳統FDTD仿真結果以及CST 仿真結果相符合,且通用HIEFDTD 方法的計算效率遠高于傳統FDTD 方法.數值計算結果表明了該程序的通用性和有效性.

2 理論與程序實現

2.1 通用HIE-FDTD 基本公式

其中,C1,C2,C3和C4為FDTD 中電場和磁場時間推進公式中的迭代系數.若令

則(1)式可以重寫為

假設目標精細結構位于z方向,則對空間求導項 ?/?z采用混合顯-隱式差分.將(2)式進行如下變形:

其中,E1=(Az1-Az2+Bz1-Bz2).忽略(3)式的末尾項,得到

采用與傳統FDTD 方法相同的Yee 元胞,并引入等效參數概念,將(5)式中的空間求導項采用二階中心差分近似,便可得到電場和磁場各個場分量的迭代公式.為避免公式表述過于冗長,這里只給出了和兩個場分量的迭代公式,其余4 個場分量的迭代公式可以參照這兩個公式給出.

構建普惠金融體系成為當前我國新一輪農村金融改革的目標,小額信貸發展是實現普惠金融的重要內容。我國小額信貸發展過程中存在的問題主要體現在對“小額信貸”的理解存在偏差,小額信用貸款公司發展速度快，不僅存在“目標偏移”的傾向，而且社會資金的商業利益與社會責任脫節的情況也時有發生。指出我國小額貸款公司發展的總體目標是實現財務績效與社會績效協調發展,為實現我國普惠金融,在路徑選擇上應該注重農村金融市場機制培育。提出從開展社會績效評價、降低交易成本等方面實現普惠金融的對策建議。

式中,Δx,Δy,Δz分別為沿x,y,z方向的空間網格大小;i,j,k為空間位置標記;n為時刻標記;Δt為時間步長.(6)式中的迭代系數分別為

其中,εx平均為等效介電系數,σx平均為等效電導率,μy平均為等效磁導系數,σmy平均為等效磁導率.等效參數的計算方法與傳統FDTD 等效參數的計算方法相同[21,22].

使用FDTD 方法模擬開域電磁問題時,需要在計算區域的截斷邊界引入吸收邊界條件.同理,使用HIE-FDTD 方法模擬開域問題,也需要引入吸收邊界.由于采用CPML 作為吸收邊界,電場和磁場分量均不需要分裂,計算公式較為簡單,因此在通用HIE-FDTD 方法中仍采用CPML 邊界條件作為吸收邊界條件.在考慮CPML 吸收邊界后,的計算公式變為

其中,

(8)式和(9)式中,Ψ為HIE-FDTD 中電場和磁場時間推進公式中卷積項的離散循環卷積形式,且CPML吸收邊界各參數的取值與其在傳統FDTD 方法中的取值一致.在非CPML 區域,kv的取值恒為1,av和bv的取值恒為0(v=x,y,z),此時(8)式退化為(6)式.引入CPML 后,電磁場其余場分量的迭代公式可參照(8)式給出[23-25].

顯然,如(8)式所示,在計算n+1 時刻的Ex時需要用到同一時刻的Hy分量.將(8b)式代入(8a)式中,整理得到

其中,

同理,在計算n+1 時刻的Ey時需要用到同一時刻的Hx分量.采用與Ex相同的處理方法,得到Ey的計算公式如下:

其中,

綜上所述,通用HIE-FDTD 方法的求解可通過如下方式完成: 首先,通過顯式迭代方程計算;再利用(10)式和(11)式,通過解三對角矩陣方程求解;在此基礎上,再通過顯式迭代方程求解.完成一次時間步迭代,通用HIE-FDTD 方法共需要求解4 個顯式迭代方程和2 個三對角矩陣方程.

2.2 程序實現

HIE-FDTD 方法與傳統FDTD 方法在程序實現時最主要的差異在于,HIE-FDTD 需要構建三對角矩陣并求解其逆矩陣.以(10)式為例,在求解時,需要沿z方向以列為單位對其進行求解,且求解每一列時均需要構建一個三對角矩陣.假設計算區域沿x方向的網格數為Nx,沿y方向的網格數為Ny,沿z方向的網格數為Nz,且Nz=7,則求解每一列所需要構建的三對角矩陣方程為

其中,τ1k和τ2k分別為k取不同值時所對應的τ1和τ2的值,Exk為沿z方向第k個待求電場量,rk為k取不同值時(10)式等號右側的已知量.顯然,完成所有求解需要構建Nx×Ny個三對角矩陣,每個三對角矩陣的維度為(Nz—1)×(Nz—1),并執行Nx×Ny次的矩陣求逆計算.當計算區域較大時,執行程序將會耗費大量的內存和時間.為解決這一問題,本文提出了一種縮減三對角矩陣數量的方法.

以圖1 所示模型為例,對通用HIE-FDTD 方法中三對角矩陣構建以及三對角矩陣數量縮減方法進行說明.如圖1 所示,模型由目標和空氣兩部分組成,其中,目標由耶路撒冷十字金屬層和長方體介質層組成.模型沿x和y方向的網格數均為18,沿z方向的網格數為10.由于本節只介紹三對角矩陣構建以及三對角矩陣數量縮減方法,因此這里只給出模型的網格尺寸,且整個模型用均勻網格線進行剖分,空間網格線如圖1 中虛線所示.

圖1 模型結構及網格分布Fig.1.Structure and grid distribution of the model.

對模型中的媒質進行編號,將空氣標記為1,金屬標記為2,介質標記為3.圖2 為模型在x=4這一截面的介質分布矩陣.

圖2 模型在x= 4 截面的介質分布矩陣Fig.2.Media distribution matrix of x= 4.

由圖2 可見,當計算x=4 這一截面的時,共需要構建18 個三對角矩陣,每個三對角矩陣的維度為9×9.

當HIE-FDTD 方法中不引入等效參數概念時,三對角矩陣的值可由沿z方向每一列的介質組合決定.相同的介質組合對應的三對角矩陣是相同的.因此可以找出模型中不同的介質組合,求出這些不同介質組合對應的三對角矩陣即可.如圖2 所示,雖然在x=4 這一截面有18 組介質組合,但多數介質組合是重復的,不重復的介質組合只有3 種,如圖2 中陰影部分所示.

當HIE-FDTD 方法中引入等效參數概念時,參照傳統FDTD 等效參數計算方法,D1,D2中的εx平均和σx平均是由(i,j,k),(i,j-1,k),(i,j,k-1)和(i,j-1,k-1) 4 個網格的ε和σ分別求平均得到,而D3,D4中的μy平均和σmy平均則是由(i,j,k)和(i,j-1,k)兩個網格的μ和σm分別求平均得到.因此,對于相同的介質組合,其對應的三對角矩陣也不一定相同.如圖2 中的y=3 和y=4兩列的介質組合相同,但其對應的三對角矩陣的值卻不同.因此,當在HIE-FDTD 方法中引入等效參數概念時,不可以再按照介質組合是否重復來縮減三對角矩陣數量.

由于三對角矩陣是基于τ1和τ2構建的,因此可以從τ1和τ2入手尋求減少三對角矩陣數量的方法.觀察τ1和τ2的計算公式可以發現,τ1和τ2的值是隨空間位置變化的.因此可以先計算出空間所有網格點的τ1和τ2,然后將其按列(沿z方向)分別存儲到二維矩陣T1和T2中.在此基礎上,令T=,則T矩陣中的每一列對應求解時所需要構建的一個三對角矩陣.T1矩陣、T2矩陣和T矩陣的結構如圖3 所示,其中,τ1-i,j,k表示空間位置(i,j,k)處的τ1值,τ2-i,j,k表示空間位置為(i,j,k)處的τ2值,Nz為計算區域沿z方向的網格數.

圖3 T1 矩陣、T2 矩陣和T 矩陣結構圖Fig.3.Structure of T1,T2 and T.

根據電磁目標的結構特征,可以預見,T矩陣中可能存在一些相同的列,基于這些相同的列所構建的三對角矩陣是完全相同的.因此可以在構建三對角矩陣之前,先剔除T矩陣中的重復列,這樣便可以有效減少三對角矩陣的數量.在每一個時間步中計算時,只需要根據空間位置標記索引到對應的三對角矩陣進行求解即可.

圖1 所示模型結構簡單,網格規模較小.實際中的模型往往結構復雜,尺寸較大[26].為說明三對角矩陣縮減方法的有效性,以圖1 模型中的目標部分為基礎,分別構建了5×5 和10×10 的目標陣列,并在目標陣列外圍各設置2 層空氣.5×5 目標陣列的網格規模為74×74×10,10×10 目標陣列的網格規模為144×144×10.分別采用直接法和縮減法對單目標模型和陣列目標模型進行三對角矩陣構建及其逆矩陣計算.兩種方法需要構建的三對角矩陣數量、程序運行時間以及程序所需內存如表1 所列.由表1 可見,采用縮減法可以有效減少三對角矩陣數量,有效降低程序運行所需內存,有效縮短程序運行時間.

表1 直接法和縮減法對比結果Table 1.Comparison results of the direct method and the reduction method.

3 數值結果

為驗證本文所提出的通用HIE-FDTD 程序實現方法的有效性,依據該方法編寫了Matlab 程序,并用程序計算了兩個數值算例,分別為平面波照射介質板和雙頻微帶倒F 天線.在這兩個例子中,目標的精細結構均位于z方向,因此對于空間求導項?/?z采用混合顯-隱式差分算法.

3.1 平面波照射介質板模型

介質板由介質基板和金屬表面涂層兩部分組成.介質基板的介質參數為εr=2.2,μr=1,σ=0.04,σm=0 .介質基板沿x和y方向的長度均為210 mm,沿z方向的厚度為0.2 mm.金屬涂層由5×5 個耶路撒冷十字組成,厚度為0.2 mm.介質板結構及耶路撒冷十字尺寸如圖4 所示.整個計算空間網格數為120×120×56.網格大小分別為:在空氣中,Δx=Δy=Δz=5 mm;在介質基板和金屬涂層中,Δx=Δy=3 mm,Δz=0.1 mm.

圖4 介質板結構及尺寸Fig.4.Structure and dimensions of the dielectric plate.

假設平面波從金屬涂層一側照射到介質板上,在介質板中軸線上設置兩個觀察點,計算觀察點處電場的時域波形.觀察點1 位于介質板有金屬涂層一側,且距金屬涂層的距離為10.2 mm.觀察點2位于介質板無金屬涂層一側,且距無金屬涂層的距離為10.2 mm.平面波的激勵源為高斯脈沖,其時域表達式為

其中,τ=3.5×10-10,t0=0.8τ.

觀察點1 處Ex的時域波形如圖5(a)所示,觀察點2 處Ex的時域波形如圖5(b)所示.從圖5 可以看到,通用HIE-FDTD 程序的計算結果與傳統FDTD 程序的計算結果相符合.從仿真時間上來看,傳統FDTD 方法的仿真時間為300 min,而通用HIE-FDTD 方法的仿真時間僅需32 min,約為傳統FDTD 方法仿真時間的1/10.此外,采用直接法的HIE-FDTD 程序的仿真時間為45 min,與通用HIE-FDTD 程序仿真時間的差異主要是因為需構建更多三對角矩陣及其逆矩陣造成的.可以預見,當模型網格規模更大時,采用直接法的HIEFDTD 程序的計算時間將會進一步增加.

圖5 觀察點處Ex 的時域波形 (a) 觀察點1;(b) 觀察點2Fig.5.Time domain waveforms of Ex at the observation point: (a) Observation point one;(b) observation point two.

3.2 雙頻微帶倒F 天線模型

雙頻微帶倒F 天線的結構如圖6 所示.金屬地板和倒F 結構分別刻蝕在介質板的兩側.介質基板的介質參數為εr=2.2,μr=1,σ=0,σm=0.介質基板沿x和y方向的長度分別為38.4 和39 mm,沿z方向的厚度為0.8 mm.整個計算空間網格數為56×54×46,且沿x和y方向的最小網格值為2.4 mm,沿z方向的最小網格值為0.2 mm.

圖6 雙頻微帶倒F 天線結構圖Fig.6.Structure of the dual-frequency microstrip inverted-F antenna.

雙頻微帶倒F 天線的S11仿真結果如圖7 所示.可以看出,倒F 天線的工作頻率為2.43 和4.97 GHz.圖7 同時給出了傳統FDTD,HIE-FDTD 和CST仿真結果,三條S11曲線符合度較好.圖8 給出了天線在2.43 GHz 的輻射方向圖,圖9 給出了天線在4.96 GHz 的輻射方向圖.從圖8 和圖9 可以看出HIE-FDTD 的仿真結果與CST 仿真結果基本一致.從仿真時間上來看,傳統FDTD 方法的仿真時間為420 s,通用HIE-FDTD 方法的仿真時間為90 s,約為傳統FDTD 方法仿真時間的1/5.此外,采用直接法的HIE-FDTD 程序的仿真時間為113 s.可見當模型網格規模較小時,縮減三對角矩陣數量方法對于總仿真時間的減少作用并不明顯.此時,仿真時間主要由時間步長和仿真步數決定.

圖7 倒F 天線S11Fig.7.S11 of the inverted-F antenna.

圖8 天線在2.43 GHz 的輻射方向圖 (a) XOY 面;(b) XOZ 面Fig.8.Radiation patterns of the inverted-F antenna at 2.43 GHz: (a) XOY plane;(b) XOZ plane.

圖9 天線在4.96 GHz 的輻射方向圖 (a) XOY 面;(b) XOZ 面Fig.9.Radiation patterns of the inverted-F antenna at 4.96 GHz: (a) XOY plane;(b) XOZ plane.

4 總結

本文提出了一種通用HIE-FDTD 方法,并依據該方法編寫了一套通用HIE-FDTD 程序.使用該套程序可以實現對線性、非色散空間中,在一維方向具有精細結構的任意電磁目標的仿真.通用HIE-FDTD 程序可實現三對角矩陣的自動構建,并對三對角矩陣數量進行自動縮減,使用一套程序即可實現對任意目標的模擬,而不是像現有HIEFDTD 程序那樣仿真一個模型便需要修改一次程序.另外,在通用HIE-FDTD 方法中,計算迭代系數時考慮了微分形式麥克斯韋方程在介質分界面失效的問題,并引入等效參數方法解決這一問題.在以往的HIE-FDTD 文獻和專著中均沒有提到等效參數方法這一概念.在通用HIE-FDTD 基本公式中,加入CPML 吸收邊界項,使得通用HIEFDTD 程序不僅可以模擬封閉區域問題,也能模擬開放區域問題.使用通用HIE-FDTD 程序仿真了平面波照射介質板模型以及雙頻微帶倒F 天線模型,仿真結果與傳統FDTD 仿真結果以及CST軟件仿真結果相符合,且通用HIE-FDTD 程序的仿真時間遠少于傳統FDTD 的仿真時間.本文提出的通用HIE-FDTD 方法只適用于非色散空間,后續研究中可結合輔助方程等方法將其擴展到色散空間,使其應用范圍更加廣泛.