汪維軍，林 歡，陳 華,2

(1.昆明理工大學(xué) 理學(xué)院，云南 昆明 650500；2.昆明理工大學(xué) 非常規(guī)冶金省部共建教育部重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650093)

20世紀(jì)60年代末，微波加熱在生物醫(yī)學(xué)工程、食品工業(yè)和橡膠工業(yè)中首先獲得成功應(yīng)用，如今，微波加熱已經(jīng)成為一種新興的加熱技術(shù)，并在冶金、材料、化學(xué)和能源等諸多領(lǐng)域受到越來越廣泛的關(guān)注[1-2]，而電磁性能的改變對(duì)物料吸收微波能力、效果也有著明顯的影響[3]；因此，研究腔體內(nèi)的電磁特性對(duì)微波加熱至關(guān)重要。1996年，Yee提出了時(shí)域有限差分法(Finite Difference Time Domain，F(xiàn)DTD)[4-5]，作為時(shí)域方法中最為經(jīng)典的算法，該方法以其原理直觀、編程簡(jiǎn)單和實(shí)用性強(qiáng)等特點(diǎn)，成為了一種典型的全波分析方法，一直在時(shí)域方法中占主導(dǎo)地位，在計(jì)算電磁學(xué)領(lǐng)域內(nèi)被廣泛、深入地研究，極大地推動(dòng)了整個(gè)時(shí)域方法的發(fā)展[6]。FDTD法直接將有限差分式代替Maxwell旋度方程中的微分式，得到關(guān)于場(chǎng)分量的有限差分式，用具有相同電參量的空間網(wǎng)格去模擬研究體，選取合適的場(chǎng)初值和計(jì)算空間的邊界條件，進(jìn)行求解[7-9]。戴國(guó)強(qiáng)等[10]通過對(duì)完美匹配層邊界條件的設(shè)置，應(yīng)用MATLAB軟件對(duì)自由空間散射場(chǎng)進(jìn)行了仿真，解決了電磁場(chǎng)散射等問題，給出了在不同時(shí)間步下的仿真結(jié)果。王立[11]基于MATLAB軟件，采用時(shí)域有限差分法對(duì)二維TM波在空間的傳播進(jìn)行了仿真，得到了電磁波在整個(gè)區(qū)域里面的動(dòng)態(tài)傳播，以便直觀、清楚地觀察電磁場(chǎng)的分布情況。本文應(yīng)用MATLAB軟件對(duì)波導(dǎo)和諧振腔內(nèi)電磁波和其能量分布進(jìn)行了模擬，同時(shí)也對(duì)此腔體內(nèi)不同介質(zhì)下的電磁場(chǎng)進(jìn)行了仿真，進(jìn)而得到了不同條件下傳輸波導(dǎo)的最好匹配頻率。

1 時(shí)域有限差分法(FDTD)的基本原理

FDTD算法將Maxwell旋度方程組進(jìn)行時(shí)間和空間的離散，得到遞推的差分方程[12]，用差分方程的解近似代替原方程的解。在求解過程中，需要保證方程組解的收斂性和穩(wěn)定性。對(duì)于無源的空間，在直角坐標(biāo)系下，當(dāng)μ、ε、σe和σm不是時(shí)間的函數(shù)時(shí)，Maxwell方程可以分解為6個(gè)標(biāo)量方程。K. S. Yee在對(duì)這2組標(biāo)量方程進(jìn)行差分離散時(shí)的巧妙之處，就在于Yee網(wǎng)格的提出及在時(shí)間和空間上同時(shí)對(duì)6個(gè)電場(chǎng)和磁場(chǎng)分量進(jìn)行差分，其結(jié)果就是使得電場(chǎng)和磁場(chǎng)在時(shí)間和空間上都相差半個(gè)網(wǎng)格。Yee網(wǎng)格如圖1所示。

圖1 FDTD離散中的Yee元胞

(1)

(2)

(3)

(4)

(5)

(6)

式中，Ex、Ey、Ez為電場(chǎng)強(qiáng)度；Hx、Hy、Hz為磁場(chǎng)強(qiáng)度；μ為媒質(zhì)的磁導(dǎo)率；ε為介電常數(shù)；σe和σm分別為媒質(zhì)的電損耗和磁損耗。

根據(jù)Yee提出的元胞原理，將式1進(jìn)行離散得到：

(7)

Δx、Δy和Δz為空間步長(zhǎng)；Δt為時(shí)間步長(zhǎng)。

2 FDTD中的激勵(lì)源

用FDTD方法分析電磁問題時(shí)一個(gè)重要的任務(wù)是對(duì)激勵(lì)源的模擬，即選擇合適的入射波形式以及用適當(dāng)?shù)姆椒▽⑷肷洳尤氲紽DTD迭代中。本文選取的是微分高斯脈沖，將高斯脈沖求導(dǎo)后得到的微分高斯脈沖函數(shù)如下：

(8)

其優(yōu)點(diǎn)是不含零頻率分量，頻域形式為：

(9)

當(dāng)以微分高斯脈沖為入射波時(shí)，如果關(guān)心的頻譜范圍的最高頻率為fmax，則對(duì)應(yīng)的真空波長(zhǎng)為λmin(λmin=c/fmax，c為自由空間波速)。若FDTD的Yee元胞尺寸取為σ=λmin/N，其中N≥10，且取時(shí)間間隔為Δt(Δt=σ/(2c))，則有：

(10)

若脈沖寬度為τ，則在1個(gè)脈沖時(shí)間內(nèi)，應(yīng)至少取采樣點(diǎn)總數(shù)Nτ為：

(11)

3 吸收邊界條件

自從Yee提出FDTD方法以來，吸收邊界條件的研究就一直沒有停止過。目前，吸收邊界條件主要有2種：一種是在邊界上引入吸收材料，電磁波在無反射地進(jìn)入吸收材料后被衰減掉，如PML(Perfectly Matched Layer)，這種方法構(gòu)造復(fù)雜，內(nèi)存需求較大，但在很大的入射角度上吸收效果較好；另一種是從外行波方程出發(fā)構(gòu)造的透射邊界條件，如Mur邊界條件等，這種類型的透射邊界條件具有構(gòu)造簡(jiǎn)單、內(nèi)存需求小和基本上不額外消耗內(nèi)存等特點(diǎn)。

由于計(jì)算機(jī)的局限性，基于有限差分的電磁數(shù)值模擬只能在有限區(qū)域內(nèi)進(jìn)行，這樣必須對(duì)計(jì)算區(qū)域的邊界進(jìn)行截?cái)嗵幚韀13]，各向異性介質(zhì)完全匹配層(Uniaxial PML，UPML)邊界條件不需要對(duì)場(chǎng)分量進(jìn)行分裂，迭代公式簡(jiǎn)單，便于程序?qū)崿F(xiàn)[14]，且UPML對(duì)各頻段電磁波的吸收具有比較好的穩(wěn)定性[15]；因此，本文中將會(huì)用到UPML吸收邊界條件、各向異性介質(zhì)麥克斯韋旋度方程(無源)在時(shí)諧場(chǎng)情形時(shí)為：

▽×H=jωε·E
▽×E=-jωμ·H

(12)

為便于得到時(shí)間推進(jìn)式，引入中間變量D、B，經(jīng)過推導(dǎo)便能得到截?cái)嘟^緣介質(zhì)的UPML三維時(shí)域方程，其時(shí)間推進(jìn)步驟為H→D→E→B→H，截?cái)嘟^緣介質(zhì)的UPML三維FDTD的推導(dǎo)及公式參見文獻(xiàn)[16]。

推導(dǎo)出FDTD的更新方程后，可以構(gòu)造出時(shí)間步進(jìn)算法。MATLAB程序流程圖如圖2所示。

圖2 FDTD算法流程圖

4 數(shù)值算例

以X波段WR90矩形波導(dǎo)為例(見圖3)，其寬為22.86 mm，高為10.16 mm，模式的理論截止頻率為6.562 GHz。應(yīng)用MATLAB軟件仿真時(shí)，使用UPML吸收邊界條件截?cái)嗑匦尾▽?dǎo)的兩端，其他的四邊用PEC截?cái)?。矩形波?dǎo)在橫截面上劃分為24×12個(gè)網(wǎng)格，而在傳播方向上大約劃分為60個(gè)網(wǎng)格。

圖3 X波段WR90矩形波導(dǎo)

應(yīng)用MATLAB軟件對(duì)該矩形波導(dǎo)的TE10模式的波阻抗進(jìn)行仿真，得到波阻抗與頻率的關(guān)系，將得到的結(jié)果與文獻(xiàn)[17]中用FDTD仿真該矩形波導(dǎo)的結(jié)果進(jìn)行比較，結(jié)果如圖4所示。從圖4中可以看出，隨著頻率的增大，波阻抗先增大后減小，當(dāng)達(dá)到最大時(shí)，截止頻率為6.56 GHz。2種仿真的結(jié)果與理論結(jié)果一樣，說明此MATLAB程序的正確性。

圖4 X波段WR90矩形波導(dǎo)的波阻抗

本文研究的模型如圖5所示。該模型由諧振腔和傳輸波導(dǎo)組成。其中，諧振腔的長(zhǎng)、寬、高尺寸分別為200、100和100 mm，傳輸波導(dǎo)為BJ—26型波導(dǎo)，其長(zhǎng)、寬、高尺寸分別為43、140和22 mm。采用的激勵(lì)源為微分高斯脈沖源，中心頻率為2.45 GHz。為建模方便，采用等間隔離散Δx=Δy=Δz=Δ=2 mm，時(shí)間步長(zhǎng)Δt=Δ/(2c)。激勵(lì)源從x-o-z面處饋入，仿真步長(zhǎng)Δt=2 000時(shí)，應(yīng)用MATLAB軟件，仿真得到真空下腔體內(nèi)部電場(chǎng)分布圖如圖6所示。

圖5 波導(dǎo)加載諧振腔結(jié)構(gòu)的FDTD設(shè)置

圖6 真空下腔體內(nèi)部電場(chǎng)分布圖

圖6a為未加入U(xiǎn)PML邊界條件，Ey在z=30σ平面內(nèi)，Δt=2 000時(shí)刻的分布圖；圖6b為加入U(xiǎn)PML邊界條件Ey在z=30σ平面內(nèi)，Δt=2 000時(shí)刻的分布圖。從圖6a可以看出，無邊界條件下電磁場(chǎng)在邊界上有反射，腔內(nèi)的電場(chǎng)比較強(qiáng)；由圖6b可知，加入U(xiǎn)PML邊界條件后，在該邊界處電磁波不再向前傳播，而是全部被吸收，腔體內(nèi)電場(chǎng)比較弱，沒有明顯的反射，表明UPML層的吸收效果較好。

上述研究了空腔下電磁場(chǎng)的分布，但是在大多數(shù)情況下，所面對(duì)的腔體里是充滿介質(zhì)的，因此，本文分別采用干燥的木頭(εr=3)和石英(εr=5)[18]對(duì)腔體進(jìn)行填充，應(yīng)用MATLAB軟件進(jìn)行仿真，得到在不同介質(zhì)、相同時(shí)間步、同一截面處電場(chǎng)的分布情況(見圖7)。

圖7 不同介質(zhì)、相同時(shí)間步、同一截面處電場(chǎng)的分布圖

從圖7可以看出，在同一時(shí)刻，在其他條件不變的情況下，隨著腔體內(nèi)所加介質(zhì)相對(duì)介電常數(shù)的增大，其腔體內(nèi)的電場(chǎng)分布越均勻，從而反映出該介質(zhì)對(duì)電磁波的吸收效果越好。在該模型加入U(xiǎn)PML邊界條件，腔體內(nèi)分別填充真空、干燥的木頭和石英介質(zhì)的條件下，可以得到y(tǒng)方向上同一截面處的s參數(shù)。當(dāng)運(yùn)行時(shí)間為2 000時(shí)間步，腔體內(nèi)分別為真空、干燥的木頭和石英介質(zhì)時(shí)，傳輸波導(dǎo)在y方向同一截面處的s11參數(shù)(即反射參量)如圖8所示。

從圖8中可以看到，s11在不同的介質(zhì)中有所不同：在空腔下匹配最好時(shí)的幅度值為-5 dB，頻率為2.42 GHz；加入干燥的木頭時(shí)的幅度值為-6.5 dB，頻率為2.35 GHz；加入石英時(shí)幅度值為-7.08 dB，頻率為2.41 GHz。說明隨著腔體內(nèi)填充介質(zhì)的相對(duì)介電常數(shù)的增大，吸收電磁波的效果增強(qiáng)。雖然在模擬中應(yīng)用了UPML吸收邊界條件，但是在仿真中存在頻差或者也有可能引入了其他的誤差，所以，在不同的介質(zhì)下達(dá)到最佳匹配時(shí)的頻率與理論上相比存在偏差。

5 結(jié)語

利用時(shí)域有限差分法，在單一激勵(lì)源的條件下，對(duì)傳輸波導(dǎo)和諧振腔內(nèi)電磁場(chǎng)的分布進(jìn)行了模擬。首先利用求矩形波導(dǎo)的波阻抗對(duì)程序進(jìn)行了驗(yàn)證；其次分別在傳輸波導(dǎo)是否加邊界條件的情況下，對(duì)波導(dǎo)內(nèi)同一平面處的電磁波做了比較；最后在腔體內(nèi)填充了不同的介質(zhì)，進(jìn)而得到了在不同的條件下腔體內(nèi)的電磁場(chǎng)的分布、電介質(zhì)對(duì)電磁波傳播的影響以及在同一截面處的散射參量。

[1] 胡兵，黃柱成，王華. 鐵及其氧化物在微波場(chǎng)中的升溫特性[J]. 中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2013, 44(8):3095-3101．

[2] 劉宏強(qiáng)，魏賢智，李飛，等. 基于射頻隱身的雷達(dá)跟蹤狀態(tài)下單次輻射能量實(shí)時(shí)控制方法[J]. 電子學(xué)報(bào)，2015,43(10):2047-2052．

[3] Yee K S. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media[J]. IEEE Trans. Antennas Propagat., 1966, 14:302-307.

[4] 王長(zhǎng)清，祝西里. 電磁場(chǎng)中的時(shí)域有限差分法[M]. 北京：北京大學(xué)出版社，1994.

[5] 杜磊. 時(shí)域有限元電磁計(jì)算方法的研究[D]. 南京：南京理工大學(xué)，2009.

[6] Jha A K, Akhtar M J. A generalized rectangular cavity approach for determination of complex permittivity of materials[J]. IEEE Transactions on Instrumentation and Measurement, 2014, 63(11):2632-2641.

[7] Guo W H, Li W J, Huang Y Z. Computation of resonant frequencies and quality factors of cavities by FDTD technique and Padé approximation[J]. IEEE Microwave and Wireless Components Letters, 2001, 11(5):223-225.

[8] Oskooi A F, Roundy D, Ibanescu M, et al. MEEP:A flexible free-software package for electromagnetic simulations by the FDTD method[J]. Computer Physics Communications, 2010, 181(3):687-702.

[9] 潘小娟，陳津，張猛，等. 微波加熱含碳碳酸錳礦粉升溫機(jī)理[J]. 中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2008, 39(6):1233-1238．

[10] 戴國(guó)強(qiáng)，余震虹. 時(shí)域有限差分方法的研究[J]. 現(xiàn)代電子技術(shù)，2013, 36(1):140-143.

[11] 王立. 時(shí)域有限差分法仿真二維電磁波傳播[J]. 艦船科學(xué)技術(shù)，2011, 33(2):49-52.

[12] Elsherbeni A, Demir V. MATLAB模擬的電磁學(xué)時(shí)域有限差分法[M]. 喻志遠(yuǎn)，譯. 北京：國(guó)防工業(yè)出版社，2013.

[13] 姜彥南，等. 瞬變電磁場(chǎng)模擬的CPML吸收邊界條件[J]. 計(jì)算物理，2015, 32(6):701-708.

[14] 魏麗君，張彬，張若冰. UPML吸收邊界條件的無截?cái)鄶?shù)值模擬[J]. 數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2016, 46(6):254-260.

[15] 于濤，戚宗鋒，李志鵬. 時(shí)域有限差分法常用吸收邊界的性能分析[J]. 海軍航空工程學(xué)院學(xué)報(bào)，2016, 31(5):506-512.

[16] 葛德彪，閆玉波. 電磁波時(shí)域有限差分方法[M]. 3版. 西安：西安電子科技大學(xué)出版社，2011.

[17] 楊儒貴. 高等電磁場(chǎng)理論[M]. 北京:高等教育出版社，2008.

[18] 謝處方，饒克謹(jǐn). 電磁場(chǎng)與電磁波[M]. 4版. 北京：高等教育出版社，2006.