基于MATLAB的微波加熱腔內介質擾動下的數值模擬研究*

汪維軍，林 歡，陳 華,2

(1.昆明理工大學 理學院，云南 昆明 650500；2.昆明理工大學 非常規冶金省部共建教育部重點實驗室，云南 昆明 650093)

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

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

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

圖1 FDTD離散中的Yee元胞

(1)

(2)

(3)

(4)

(5)

(6)

式中，Ex、Ey、Ez為電場強度；Hx、Hy、Hz為磁場強度；μ為媒質的磁導率；ε為介電常數；σe和σm分別為媒質的電損耗和磁損耗。

根據Yee提出的元胞原理，將式1進行離散得到：

(7)

Δx、Δy和Δz為空間步長；Δt為時間步長。

2 FDTD中的激勵源

用FDTD方法分析電磁問題時一個重要的任務是對激勵源的模擬，即選擇合適的入射波形式以及用適當的方法將入射波加入到FDTD迭代中。本文選取的是微分高斯脈沖，將高斯脈沖求導后得到的微分高斯脈沖函數如下：

(8)

其優點是不含零頻率分量，頻域形式為：

(9)

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

(10)

若脈沖寬度為τ，則在1個脈沖時間內，應至少取采樣點總數Nτ為：

(11)

3 吸收邊界條件

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

由于計算機的局限性，基于有限差分的電磁數值模擬只能在有限區域內進行，這樣必須對計算區域的邊界進行截斷處理[13]，各向異性介質完全匹配層(Uniaxial PML，UPML)邊界條件不需要對場分量進行分裂，迭代公式簡單，便于程序實現[14]，且UPML對各頻段電磁波的吸收具有比較好的穩定性[15]；因此，本文中將會用到UPML吸收邊界條件、各向異性介質麥克斯韋旋度方程(無源)在時諧場情形時為：

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

(12)

為便于得到時間推進式，引入中間變量D、B，經過推導便能得到截斷絕緣介質的UPML三維時域方程，其時間推進步驟為H→D→E→B→H，截斷絕緣介質的UPML三維FDTD的推導及公式參見文獻[16]。

推導出FDTD的更新方程后，可以構造出時間步進算法。MATLAB程序流程圖如圖2所示。

圖2 FDTD算法流程圖

4 數值算例

以X波段WR90矩形波導為例(見圖3)，其寬為22.86 mm，高為10.16 mm，模式的理論截止頻率為6.562 GHz。應用MATLAB軟件仿真時，使用UPML吸收邊界條件截斷矩形波導的兩端，其他的四邊用PEC截斷。矩形波導在橫截面上劃分為24×12個網格，而在傳播方向上大約劃分為60個網格。

圖3 X波段WR90矩形波導

應用MATLAB軟件對該矩形波導的TE10模式的波阻抗進行仿真，得到波阻抗與頻率的關系，將得到的結果與文獻[17]中用FDTD仿真該矩形波導的結果進行比較，結果如圖4所示。從圖4中可以看出，隨著頻率的增大，波阻抗先增大后減小，當達到最大時，截止頻率為6.56 GHz。2種仿真的結果與理論結果一樣，說明此MATLAB程序的正確性。

圖4 X波段WR90矩形波導的波阻抗

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

圖5 波導加載諧振腔結構的FDTD設置

圖6 真空下腔體內部電場分布圖

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

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

圖7 不同介質、相同時間步、同一截面處電場的分布圖

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

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

5 結語

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

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

[2] 劉宏強，魏賢智，李飛，等. 基于射頻隱身的雷達跟蹤狀態下單次輻射能量實時控制方法[J]. 電子學報，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] 王長清，祝西里. 電磁場中的時域有限差分法[M]. 北京：北京大學出版社，1994.

[5] 杜磊. 時域有限元電磁計算方法的研究[D]. 南京：南京理工大學，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]. 中南大學學報：自然科學版，2008, 39(6):1233-1238．

[10] 戴國強，余震虹. 時域有限差分方法的研究[J]. 現代電子技術，2013, 36(1):140-143.

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

[12] Elsherbeni A, Demir V. MATLAB模擬的電磁學時域有限差分法[M]. 喻志遠，譯. 北京：國防工業出版社，2013.

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

[14] 魏麗君，張彬，張若冰. UPML吸收邊界條件的無截斷數值模擬[J]. 數學的實踐與認識，2016, 46(6):254-260.

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

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

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

[18] 謝處方，饒克謹. 電磁場與電磁波[M]. 4版. 北京：高等教育出版社，2006.