聶睿瑞 李 勃

(1.南京信息職業(yè)技術學院智能交通學院，江蘇 南京 210023；2.南京航空航天大學無人機研究院，江蘇 南京 210016)

在現(xiàn)有的電磁場算法中，時域有限差分法FDTD的使用越來越廣泛[1－4]。因為沒有引入數(shù)學計算，傳統(tǒng)FDTD 運算法則的網(wǎng)格劃分較為簡單。然而，需要解決曲面問題時，傳統(tǒng)FDTD 算法就不太適用了。對于電氣設備的不連續(xù)曲線表面和多重材料結構，傳統(tǒng)FDTD 算法還會產(chǎn)生發(fā)散和不穩(wěn)定現(xiàn)象。所以，很多項目都把FDTD 算法的穩(wěn)定性作為研究重點[5－7]。

本文提出了一種基于非標準曲線算子的三維高階FDTD 算法，可以用于解決復雜介質結構電氣設備的電磁兼容問題。該三維高階FDTD 算法使用了足夠多的點來表示空間微商，而不是像傳統(tǒng)FDTD算法那樣只取兩個點。在完全導體邊界或吸收邊界處，空間模板不可避免的會擴大，采用緊湊的曲線差分方案，能將無限空間轉換為修正過的PML 吸收邊界條件[8－9]?？紤]到那些界面和任何網(wǎng)格都不相符的材料，該算法使用了特殊的轉換法。通過對多種三維曲線以及不同介質結構的電磁兼容性計算，結果表明該方法是有效的。

1 高階非標準FDTD 算法

1.1 非標準形狀的結構

三維高階FDTD 算法采用式(1)來描述三維電磁場的結構[10]:

式中:g是三維坐標函數(shù)；CS是修正函數(shù)；Rm和Pm是穩(wěn)定性參數(shù)；Dq是最佳節(jié)點算子。q是常規(guī)坐標系(u，v，w)的變量，M表示準確度的階數(shù)。該差分法能有效降低色差的影響。而且，只要參數(shù)和滿足式(2)和式(3)，就可以確定幾何細節(jié)及場量分配的穩(wěn)定性。

式中:L因子表示模板的數(shù)量，一般取L＝3，lδq表示各軸上的坐標。修正函數(shù)CS(klδq)用于確保從連續(xù)到離散的平滑過渡，其幅度取決于波數(shù)k，對網(wǎng)格電磁向量進行傅立葉變換可以得出，CS一般由式(4)求得:

式中:rA、rB是三維計算因子。

同樣可以推出在方向u和方向w上的方程。

高階非標準算子可以應用到麥克斯韋旋度方程中的空間微商和時間微商的離散化中。根據(jù)其近似程度，通過一組參數(shù)關系式能夠得出基本模型。

空間微商:

時間微商:

式中:?m是微分參數(shù)，CT(δt)是T[·]的修正函數(shù)，b1和b2是特定的調節(jié)參數(shù)。除了全導電墻，式(7)在域中處處適用。對于PEC wall，其模板在格點的兩側都會至少延伸兩個節(jié)點。相應的曲線壓縮公式為:

式中:τA，τB，τC是三維離散積分系數(shù)；S是將曲線從三維壓縮到二維的參數(shù)。

此方程可自適應于曲率的變化，對連續(xù)物理空間進行離散化，轉換為二重拓撲網(wǎng)格結構。

1.2 非正交二重網(wǎng)格

將給定三維區(qū)域按照坐標系(u＝iδu，v＝jδv，w＝kδw)進行劃分，如圖1 所示。各個主單元的中心是(i，j，k)，而第二重單元的中心則是主單元的頂點。磁場向量H(位于主表面的中心)的共變量hq以及逆變量hp和電場向量E(位于邊緣的中心)的共變量eq以及逆變量ep相互交錯。

圖1 FDTD 二重網(wǎng)格劃分圖

我們是根據(jù)表面的場通量可定義為f(p)＝g1/2f(p)(f是電場或磁場)來進行這樣的劃分，其中gpq為坐標系的度量。由FDTD 原理可知，fq＝gqp fp且fp＝gpqfq。再引入線性算子Q(p)，得到f(p)＝Q(p)[fq]，式中用到原有分量fq及相鄰分量fq＋1和fq＋2(q＋1、q＋2表示u、v、w的連續(xù)循環(huán))。比如u軸上Q(u)[hu]的表達式是:

式中:Θpq＝g1/2gpq，因為省略了二重交叉場復雜的投影分量，這樣對通量的分析就相對簡便。通過以上步驟并替換旋度算子，可將安培和法拉第定律寫為:

上式中Ecv＝[eu ev ew]和Hcv＝[hu hv hw]是未知共變分量的矩陣；L＝[Lu，Lv，Lw]為空間算子；J＝σE和M＝σ?H是對應的傳導電流密度；Yt和Rt是定義每個介質表面的基本矩陣；GH，GE是合適的度量張量，μ是磁導率。矩陣TE和TH根據(jù)T[·]集合了所有的高階臨時微商，且容許修正。因為式(10)和(11)把所有的高階臨時微商都聚集到一個矩陣中，這樣就能利用臨時儲存向量通過修正廣義跳點法把每個時間步長分成和階數(shù)相同的數(shù)量。

采用馮諾依曼分析法來分析式(6)和式(7)的穩(wěn)定性，可知:

可調節(jié)的準確度系數(shù)M和L極大地改善了色散關系。把高階非標準FDTD 算法的色散關系和二階Yee 單元進行比較，經(jīng)數(shù)學運算有:

由此看出上式只是普通FDTD 關系FFDTD(·)中很小的一部分。例如，取M＝4，L＝3，c′是相速度的數(shù)值，β＝2π/(kδw)，θ為入射角時，有:

式中:nst的含義是非標準。

合適的和分量取值使得色散關系得到了明顯提高。高階非標準FDTD 算法能夠改善后續(xù)計算的不穩(wěn)定性、減小色散誤差且無需使用共形技術?；谡w場的實現(xiàn)，發(fā)散的問題可采用曲線非標準PML 邊界條件來解決。選用合理的縮放比例，保持原場變化，進行合適的坐標變換，就能達到最優(yōu)吸收邊界條件。

1.3 空間低通濾波過程

由于對介質表面采取了不定的調整，高階非標準FDTD 算法的不穩(wěn)定因素就是不均勻的網(wǎng)格劃分方式。而消除寄生高頻成分是一個可行的應對辦法。因此，我們可以對單元平均值采取后期處理:

式中:單元邊界[u1，u2]、[v1，v2]和[w1，w2]的維數(shù)表示為[q1，q2]集合，G是向量Ecv和Hcv的任意分量。對于多重空間，該式可以在網(wǎng)格的u、v、w方向上逐一應用。再引入自由度的附加度Ξ來控制濾波，則域的內(nèi)部可記為:

式(17)的頻率響應為:

由式(17)的均衡性可知，F(xiàn)S(k)為實數(shù)，濾波只改變場的幅度。為了能夠抑制最高頻，要求FS(k)＝0。一般取Ξ的值為4，且ζf＝0.351 4、a1＝1.287 587 493 5、a2＝－1.923 789 646 3、a3＝2.575 179 324 9、a4＝0.346 821 379 6。

1.4 介質跳變邊界的收斂處理

如果對曲線介質表面的不連續(xù)性沒有采用合適的跳變條件，就會產(chǎn)生不同的相速度和波長，這會明顯影響FDTD 算法的穩(wěn)定性和收斂性。并且，電磁場在這些邊界上產(chǎn)生不連續(xù)，即便對該區(qū)域采用平滑處理方法，最終的模擬結果仍然會出現(xiàn)極大的偏差。

對于圖2 中所示的介質曲面，?。?，，)T為單位向量。任意兩個區(qū)域(mat ＝A，B)中，如果決定了分辨度就定義了特定的空間值。由此，Emat和Hmat向量就能通過合適的切線或普通連續(xù)條件從該介質的εmat和μmat中得出。選用一系列的協(xié)變分量，表面?w、hu的修正為:

圖2 介質曲面處理

式中:hu可以通過下式算出:

式(20)的三種普通變形推算為:

對于e分量而言，表達式也類似。對式(22)、(23)的非變量有:

式(19)～(25)預先修正了高階非標準FDTD 算法，這使麥克斯韋方程在整個域中都是可解的，從而大大提高了對曲線表面的計算準確度。

2 試驗與仿真

設計圖3 中所示的某機電設備電路板，將集成電路安裝在第一層板上，它們由相同的介質構成。介質參數(shù)分別為:εA＝2.5ε0，μA＝1.76μ0，σA＝0.22 S/m?；膮?shù)是:εB＝4.7ε0，μB＝μ0，σB＝0.03 S/m。對于該機電設備電路板，在建模時需要考慮損耗以及0.01 F 的去耦電容在回路中所產(chǎn)生的高頻能量。該電路板尺寸為:b1＝18.65 mm，b2＝16.93 mm，b3＝2.5 mm。把整個區(qū)域劃分為36×24×8 個單元，取δx＝δy＝δz＝1.5 mm，δt＝25.32 ps。

圖3 試驗用某機電設備電路板

圖4 和圖5 分別顯示出參量S11(端口1 上的輻射)和S21(兩個端口之間的輻射)隨頻率變化時，二階FDTD 算法和高階非標準FDTD 算法的仿真結果。從圖4 中可以看出，高階非標準FDTD 算法的準確度非常高，并且沒有出現(xiàn)任何色散情況。相反的，盡管Yee 算法采用了標準的82×64×156 網(wǎng)格，卻不能計算出準確的峰值。

圖4 S11的幅度

圖5 S21的幅度

最后，圖6 中顯示出不同的離散化程度下，在電路板某固定位置處的標準化相速度。

圖6 多種FDTD 方法計算的相速度對比

3 結論

結果表明，本文所提出的基于非標準曲線算子的三維高階FDTD 算法，其計算結果比二階FDTD算法要精確。因為隨著分辨度的下降，二階FDTD會嚴重偏離實際值。所以在分析較大結構的EMC問題時，采用高階FDTD 算法能夠較好地彌補傳統(tǒng)方法在準確度方面的缺陷。此外，高階非標準FDTD 算法還可以節(jié)省大約85%的CPU 及內(nèi)存。