聶睿瑞 李 勃

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

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

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

1 高階非標(biāo)準(zhǔn)FDTD 算法

1.1 非標(biāo)準(zhǔn)形狀的結(jié)構(gòu)

三維高階FDTD 算法采用式(1)來(lái)描述三維電磁場(chǎng)的結(jié)構(gòu)[10]:

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

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

式中:rA、rB是三維計(jì)算因子。

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

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

空間微商:

時(shí)間微商:

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

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

此方程可自適應(yīng)于曲率的變化，對(duì)連續(xù)物理空間進(jìn)行離散化，轉(zhuǎn)換為二重拓?fù)渚W(wǎng)格結(jié)構(gòu)。

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

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

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

我們是根據(jù)表面的場(chǎng)通量可定義為f(p)＝g1/2f(p)(f是電場(chǎng)或磁場(chǎng))來(lái)進(jìn)行這樣的劃分，其中g(shù)pq為坐標(biāo)系的度量。由FDTD 原理可知，fq＝gqp fp且fp＝gpqfq。再引入線(xiàn)性算子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]的表達(dá)式是:

式中:Θpq＝g1/2gpq，因?yàn)槭÷粤硕亟徊鎴?chǎng)復(fù)雜的投影分量，這樣對(duì)通量的分析就相對(duì)簡(jiǎn)便。通過(guò)以上步驟并替換旋度算子，可將安培和法拉第定律寫(xiě)為:

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

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

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

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

式中:nst的含義是非標(biāo)準(zhǔn)。

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

1.3 空間低通濾波過(guò)程

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

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

式(17)的頻率響應(yīng)為:

由式(17)的均衡性可知，F(xiàn)S(k)為實(shí)數(shù)，濾波只改變場(chǎng)的幅度。為了能夠抑制最高頻，要求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 介質(zhì)跳變邊界的收斂處理

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

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

圖2 介質(zhì)曲面處理

式中:hu可以通過(guò)下式算出:

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

對(duì)于e分量而言，表達(dá)式也類(lèi)似。對(duì)式(22)、(23)的非變量有:

式(19)～(25)預(yù)先修正了高階非標(biāo)準(zhǔn)FDTD 算法，這使麥克斯韋方程在整個(gè)域中都是可解的，從而大大提高了對(duì)曲線(xiàn)表面的計(jì)算準(zhǔn)確度。

2 試驗(yàn)與仿真

設(shè)計(jì)圖3 中所示的某機(jī)電設(shè)備電路板，將集成電路安裝在第一層板上，它們由相同的介質(zhì)構(gòu)成。介質(zhì)參數(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。對(duì)于該機(jī)電設(shè)備電路板，在建模時(shí)需要考慮損耗以及0.01 F 的去耦電容在回路中所產(chǎn)生的高頻能量。該電路板尺寸為:b1＝18.65 mm，b2＝16.93 mm，b3＝2.5 mm。把整個(gè)區(qū)域劃分為36×24×8 個(gè)單元，取δx＝δy＝δz＝1.5 mm，δt＝25.32 ps。

圖3 試驗(yàn)用某機(jī)電設(shè)備電路板

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

圖4 S11的幅度

圖5 S21的幅度

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

圖6 多種FDTD 方法計(jì)算的相速度對(duì)比

3 結(jié)論

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