有限厚度平面周期結構的快速算法研究

車永星 侯新宇

(電子薄膜與集成器件國家重點實驗室,四川 成都610054)

1.引 言

平面周期結構,尤指頻率選擇表面(Frequency Selective Surface,FSS)是由無源諧振單元按一定的排列方式組成的準二維平面周期陣列[1]。它對某些頻帶內(nèi)的入射電磁波呈全反射特性,而對另一些頻帶內(nèi)的入射波則呈透射特性。其在微波毫米波領域的應用非常廣泛,如濾波器設計、微帶電路、天線散射等。

在傳統(tǒng)的許多含有縫隙的FSS輻射結構分析時一般都假設金屬屏的厚度為零,由于新的應用以及頻率向更高發(fā)展的需要,必須將金屬FSS屏的厚度考慮在內(nèi)[2]。有許多全波法可以求解一般的孔徑問題,但缺乏簡便性和可用于零厚度孔徑下直接積分公式模型(即2.5-D模型)的內(nèi)部一致性。因此,研究可精確求解有限厚度金屬屏,且同時具有與零厚度情況下相似的簡便性的有效近似法是非常有必要的。

根據(jù)已有的關于求解周期結構問題的文獻[3][4],可以看到基于矩量法求解的積分公式法(Integral Equation-Method of Moments,IE-MoM)是一種高精度的算法。但IE-MoM算法需要可靠準確的相關格林函數(shù),而對于周期問題精確定義的格林函數(shù)會引出另一問題,即周期格林函數(shù)的慢收斂性。這就要求對現(xiàn)有的一般或特殊形式的加速技術作進一步的研究。

由文獻[5][6]可知非線性級數(shù)加速法Shanks變換在通常情況下有較好的收斂,但在靠近源的區(qū)域會遇到困難。由文獻[7]可知特殊加速算法Ewald變換則對自由空間格林函數(shù)有非?？斓氖諗?。

論文綜合采用了一般分析方法中Shanks變換和特殊分析方法中Ewald變換的優(yōu)點,針對加速多層周期格林函數(shù),尤其是矢量位格林函數(shù),給出了一種新的有效的加速算法。此算法將屏的有限厚度并入零厚度情況下的算法或公式中,同時不增加計算時間和計算復雜度。最后,文中利用此方法對有限厚度的四腿加載單元FSS做了相應計算,給出了相應的計算結果。

2.有限厚度FSS分析

2.1 積分方程-矩量法(IE-MoM)

首先,對FSS做離散化處理。在FSS上定義一個周期單元,并將每個周期單元表面離散為多個網(wǎng)格,為了較好地近似實際表面,采用三角形或矩形組成的劃分網(wǎng)格。任何由周期結構的散射或輻射所引出的問題均可以表達為如下積分方程

式中,算子L為積分-差分算子,其將未知的表面電流 a和已知的激勵v聯(lián)系在一起。在矩量法(MoM)中,將未知電流a展開為一組N維矢量基函數(shù)fi的線性組合,其表示如下

然后將積分方程(1)與矢量權函數(shù)組wi做內(nèi)積,方程可離散化為

方程(3)也可以簡單表達如下

式中 ：Z=[zij]N×N是 MoM 矩陣 ;V=[vj]N×1為激勵矢量;I=[ai]N×1表示未知電流系數(shù)。在這里權函數(shù)wi=f i.

由于電流和磁流結合不適于模擬周期結構,因此,在這里以位格林函數(shù)作為基函數(shù)和權函數(shù)。圖1為單層介質襯底上的平面周期結構,其中介質襯底厚度為h,任意形狀的相同單元分別沿x軸和y軸以周期a、b分布,平面波以角度 θin和 φin照射,在這里假設導體厚度可忽略不計。

圖1 單層介質襯底的平面FSS

如果假設表示無周期邊界的二維自由空間或層介質中的格林函數(shù)為G0(x-x′,y-y′),則對于上述周期結構的格林函數(shù)的“空間域”公式表示如下

式中：Gmn(x,y)=G0(x-x′-ma,y-y′-nb)=ksinθicosφi和=ksinθisinφi分別為波傳播矢量的x軸分量和y軸分量;m,n分別表示沿x軸方向的第m個單元,沿y軸方向的第n個單元。

應用二維泊松求和公式[7]及傅里葉變換,由公式(5)即可得到頻域的周期格林函數(shù)如下

2.2 加速算法

在矩量法求解積分方程方法中,要求有可靠穩(wěn)定的周期格林函數(shù),但這同時會引出另一問題,即周期格林函數(shù)的慢收斂。為了保證計算時間在可接受范圍內(nèi),并可以提供精確的計算結果,在計算中必須使用加速算法。一種針對加速多層周期矢量位格林函數(shù)的加速算法如下。首先將周期格林函數(shù)表達式(6)重新寫為

式中：矢量位周期格林函數(shù)由兩部分組成,其中前一項為“動態(tài)”部分,后一項為“靜態(tài)”部分。在模式和的“動態(tài)”部分的收斂中應用改進的Shanks變換

初始條件為

式中：Sn為級數(shù)的n項部分和,當n→∞時,Sn→S,S為級數(shù)和。利用上述級數(shù)收斂變換,模式和的“動態(tài)”部分將只需計算偶數(shù)階部分,且有更好的數(shù)值穩(wěn)定性和精度,得到快速收斂。在“靜態(tài)”部分的收斂問題中應用Ewald變換可得

式中的補余誤差函數(shù)將使得級數(shù)迅速收斂。將式(8)和式(10)合并后,即可實現(xiàn)對矢量位周期格林函數(shù)的加速收斂。

2.3 有限厚度FSS的算法實現(xiàn)

對于有限厚度的縫隙型FSS,可以分解為兩部分分別求解,即內(nèi)問題(縫隙內(nèi)部區(qū)域)和外問題。外問題應用2.1節(jié)給出IE-MoM 方法求解,并利用2.2節(jié)給出的加速算法加速MoM 的計算速度。對于內(nèi)問題,將縫隙腔體由具有相同幾何截面的2-D目標代替[8-10],然后通過修改相應格林函數(shù)引入縫隙厚度,利用文獻[11]給出的從有限厚度到零厚度位格林函數(shù)的平滑轉換,即可在不增加額外復雜度的情況下對縫隙內(nèi)部區(qū)域進行求解。

3.數(shù)值結果與分析

3.1 模 型

以圖2所示的四腿加載單元FSS為例。其中單元周期D=7 mm,單元縫隙寬度w=0.2 mm,偶極子長度L=4.2 mm,偶極子寬度t=1.2 mm,FSS屏厚度為h,襯底材料相對介電常數(shù)εr=4.2,損耗正切 tanδ=0.02,襯底厚度h1=2 mm.

圖2 FSS單元模型

3.2 傳輸系數(shù)

分別利用上述算法和全波電磁仿真軟件HFSS對圖2所給FSS單元的傳輸系數(shù)進行計算。假定平面波沿法向入射,在FSS屏厚度為h=0 mm,0.2 mm(大約為 λ/200,由文獻[1]可知,當屏厚度大于λ/1000時不可忽略)時的傳輸曲線如圖3所示。由圖可知,由于厚度的存在FSS屏的諧振點略向高頻偏移。通過對用文中所給方法和HFSS仿真結果的比較,兩種方法取得了很好的一致性。

3.3 遠場散射

為了檢驗厚度對遠場參數(shù)的影響,對圖2所示FSS結構的散射進行了計算。無限周期格林函數(shù)作為3×3單元陣列的一次近似,使用平面波源沿法向入射。格林函數(shù)的無限求和在前三個部分求和后截斷。輻射方向圖如圖4所示。實際上,無限方法對于無限陣列中的單元能給出更好近似解,而在有限陣列中會受到邊緣散射的影響。從圖中可以明顯地看出,本文方法與全波方法計算所得結果基本一致。

3.4 表面電場分布

圖5給出了由IE-MoM算法結合矢量位周期格林函數(shù)加速算法計算得到的FSS單元電場二維分布圖。屏厚度為h=0.2 mm,由電場方向沿y軸的平面波激勵。在兩個方向上,周期格林函數(shù)均只取前40項,計算速度得到明顯提高。

3.5 計算時間

由全波方法、一般IE-MoM算法和結合了加速算法的IE-MoM算法所得計算時間隨未知量的變化曲線如圖6所示。平面波沿法向入射,電場沿y軸方向極化。由圖可知,隨著未知量個數(shù)的線性增長,全波算法和一般IE-MoM算法所耗時間的增長率逐漸增大,而加速計算的IE-MoM算法所耗時間只近似于線性增長,雖然計算時間沒有得到大幅度的提高,但仍較另兩種方法更簡單快速。

4.結 論

將積分方程矩量法應用于有限厚度平面周期結構的電磁特性的求解計算,并綜合采用了一般加速方法Shanks變換和特殊加速方法Ewald變換的優(yōu)點,給出了針對矢量位周期格林函數(shù)加速的新的有效加速算法。通過對有限厚度屏的反射系數(shù)、傳輸系數(shù)、遠場輻射、表面電場分布的計算,表明文中所述算法與全波電磁仿真軟件 HFSS所得結果具有很好的一致性,同時,在計算時效上得到了顯著提高。

[1] M UNK B.A.頻率選擇表面理論與設計[M].侯新宇,譯.北京：科學出版社,2009.

[2] PARK H H,EOM H J.Electromagnetic scattering from multiple rectangular apertures in a thick conducting screen[J].IEEE Trans.Antennas Propagation,1999,47(6)：1056-1060.

[3] ITOH T.Numerical Techniques for Microwave and Millimeter-wave Passive Structures[M?.New York：Wiley,1989.

[4] 楊 卓,薛正輝,李偉明.平面波展開結合矩量法分析頻率選擇表面[J].電波科學學報,2009,24(2)：314-317.YANG Zhuo,XUE Zhenghui,LI Weiming.Analyzing frequency selective surface with plane wave approach and the method of moments[J].Chinese Journal of Radio Science,2009,24(2)：314-317.(in Chinese)

[5] BREZINSKI C,ZAGLIA M R.Extrapolation Methods-theory and Practice[M].Amsterdam：Elsevier Science Publishers,1991.

[6] KINAYMAN N,.AKSUN M I.Comparativestudy of acceleration techniques for integrals and series in electromagnetic problems[C]//Antennas and Propagation Society International Symposium,1995,2(2)：1037-1040.

[7] PAPOULIS A.Systems and Transforms with Applications in Optics[M].Malabar,Florida：Kreiger,1981.

[8] JORDAN K,RICHTER G,SHENG P.An efficient numerical evaluation of the Green's function for the Helmholtz operator on periodic structures[J].Journal of Computational Physics,1986,63：222-235.

[9] STEVANOVI C I,CRESPO-VALERO P,MOSIG J R.An integral equation technique for solving thick irises in rectangular waveguides[J].IEEE Trans.Microwave Theory Tech.,2006,54(1)：189-197.

[10] STEVANOVI C I,MOSIG J R.Efficient electromagnetic analysis of line-fed aperture antennas thick conducting screens[J].IEEE Trans.Antennas Propagation,2004,52(11)：2896-2903.

[11] MOSIG J R.Scattering by arbitrarily-shaped slots on thick metallic screens：An approximate solution[J].IEEE Trans.Antennas Propagation,2004,52(8)：2109-2117.