時域有限差分法連接邊界電磁泄漏

姬金祖

(北京航空航天大學 航空科學與工程學院，北京100191)

黃大慶

(中國航空工業集團公司北京航空材料研究院，北京100095)

黃沛霖 魯振毅

(北京航空航天大學 航空科學與工程學院，北京100191)

飛行器隱身技術是現代戰爭中的一項重要軍事技術，對于提高飛行器的戰場生存力意義重大［1－2］.時域有限差分法(FDTD，Finite－Difference Time－Domain)是計算金屬和復雜介質目標電磁散射、傳播特性的一種重要方法，在隱身飛機概念設計、吸波材料設計等方面具有重要作用［3－5］.

用FDTD計算目標電磁散射的過程中，一般要把計算區域分成總場區和散射場區兩個區域［6］.總場區的場量包括入射場和散射場，散射場區的場量只有散射場，沒有入射場，因此由散射場區進行遠場外推即可得到目標的雷達散射截面RCS(Radar Cross Section)［7］.總場區和散射場區的邊界為連接邊界，入射場通過連接邊界引入到總場區.

散射體位于總場區內，當總場區不存在散射目標時，理想情況下散射場區的場量為0，總場區只有入射場.但在實際計算過程中，即使總場區沒有散射目標，散射場區的場一般也不會為0，這是因為在連接邊界引入入射場時的數值計算產生了誤差，導致電磁泄漏.

文獻［8］采用分裂平面波時域有限差分法(SP－FDTD，Splitting Plane Wave FDTD)引入入射場，在一維網格中將場量分裂，建立一維輔助網格和三維網格的色散關系.文獻［9］在電磁散射計算中使用交錯布置的方法引入入射波，該方法基于分裂場麥克斯韋方程組，擴展了多點輔助入射場的作用.文獻［10］采用一維輔助多點模型，通過CN(Crank－Nicholson)差分方法引入入射波.該方法將二維色散方程投影到一維，使電磁泄漏在時間步長取CFL(Courant Friedrichs Levy)時達到很小.

實際計算表明，時間步長以及電磁波入射方向都對電磁泄漏有一定影響，本文研究這些參數對電磁泄漏的影響的大小.

1 FDTD計算區域劃分

計算目標電磁散射過程中，FDTD計算區域一般劃分為總場區和散射場區，在散射場區還設置輻射邊界，用于將近場外推至遠場，進一步計算RCS.FDTD計算區域劃分如圖1所示.

計算過程中，在連接邊界處引入入射場，在此邊界上的時間迭代公式要做一些修改，以保證該點周圍的環量均為散射場或均為總場.時間迭代結束后，散射場區只有散射場.在輻射邊界根據等效原理，將電磁場外推至遠場，得到目標的RCS.

圖1 FDTD計算區域劃分

當總場區沒有目標時，散射場區仍然會有電磁場，這是數值計算誤差引入的電磁泄漏.本文對一維和二維情形采取不同的研究方法.一維情形下，用散射場區電場的平方和來衡量電磁泄漏的大小.在二維情形下，將散射場區的電磁場通過遠場外推得到“無目標”下的遠場RCS，用這種方法得到的RCS類似于微波暗室的背景電平.研究不同時間步長、入射角度下的遠場RCS來研究電磁泄漏的程度.

2 一維情形電磁泄漏

2.1 差分方程

一維情形下電磁場微分方程具有最簡單的形式，設電磁波傳播方向為z，電場方向為x，磁場方向為y，則電磁波在真空中傳播的微分方程如式(1)所示.

式中，Ex為電場強度;Hy為磁場強度;μ0為真空磁導率;ε0為真空介電常數.設c為真空中光速，令，結合真空中光速公式，將上式改寫成較簡單的形式，如式(2)所示.

式(2)顯得更加簡單、對稱，本文涉及的電場、磁場也均采用上式簡單、對稱的形式.

用一階Mur吸收邊界條件在兩個端點進行截斷.采用一階差分格式，對微分方程進行離散.將計算區域劃分成100個網格，在距邊界20個網格處設置連接邊界條件.將ex在網格離散點處采樣，hy在網格中點處采樣，則得到101個ex離散點和100個 hy離散點，設離散點為 ex(i，n)和hy(i，n)，其中，i表示網格點編號，n表示時間步.設入射波為正諧波，波長為λ.

設網格長度為Δz，時間步長為Δt，定義空間網格因子為 fs=z/λ，時間網格因子為 ft=Δz/(cΔt).根據色散條件和穩定性條件的要求，fs應大于15，一維情形下ft應大于1.fs主要與色散誤差相關，而本文主要研究時間網格因子ft的影響，故fs一律取為20.

用fs、ft表示真空中電磁波傳播的差分格式，如式(3)所示.

該差分格式將介電常數、磁導率、網格尺寸、時間步長等參數用空間網格因子和時間網格因子進行歸一化，形式更加簡潔，便于編寫程序和進行相應的計算.

2.2 時間網格因子影響

設電磁波傳播方向為從左向右，設置不同時間網格因子，對比計算結果.運行1000個時間步后，兩種情形下電場分布對比如圖2所示.

圖2 不同時間網格因子電場分布

由圖2可見，ft的取值對電磁泄漏影響嚴重.ft=1時，電磁泄漏最小，散射場區電場強度基本為0.而ft=10時，電磁泄漏導致散射場區電場強度有較大起伏.

為定量研究ft對電磁泄漏的影響，以散射場區的電場平方和來表征泄漏的能量.取ft的值從1到10，間隔為1，分別計算其電磁泄漏能量.電磁泄漏能量隨時間網格因子的變化如圖3所示.

由計算結果可得，當 ft=1時，泄漏能量為 1.5179×10－29，非常接近 0，可以認為此時泄漏能量很小，可以忽略不計.由圖3可得，泄漏能量基本隨ft的增加而增加.當ft＜4時，泄漏能量增加很快，但當ft＞4時，能量泄漏增加很慢，基本保持在0.06 到0.07 之間.

圖3 電磁泄漏能量隨時間網格因子的變化

下面研究ft在1到2之間變化時，電磁泄漏的變化情況.ft步長取為0.1，計算結果如圖4所示.

圖4 電磁泄漏隨時間網格因子的變化

由圖4可見，時間網格因子在1到2之間變化時，電磁泄漏能量基本上線性增加.

可見，在一維FDTD中，時間網格因子在滿足穩定性條件的基礎上，取值應盡量小，即接近于1.

3 二維情形電磁泄漏

3.1 二維情形電磁泄漏比較

下面研究二維情形下的電磁泄漏.為保證穩定性，二維情形ft的下限是計算區域網格數為160×160，連接邊界到吸收邊界距離為20個網格，外推邊界到吸收邊界為10個網格.fs與一維情形相同，取為20.電磁波在真空中傳播，無散射體時，不涉及到邊界條件，故橫磁波 TM(Transverse Magnetic)和橫電波TE(Transverse Electric)沒有本質區別.本文以TM波為例來進行分析研究.

入射角度定義為入射波來波方向與水平線夾角，逆時針為正，如圖5所示.

ft取1.5 和10 兩個值，φ 取0°和45°，計算不存在散射目標時的遠場雙站RCS，得到4條曲線.雙站角為0°～359°，步長1°.將4條曲線在一張圖上進行對比，如圖6所示.二維RCS單位為m，在圖中換算為對數值，單位為dBm

圖5 入射角度示意圖

圖6 無目標下遠場雙站RCS

上面計算的RCS可以理解為類似微波暗室的“背景電平”.由上圖可見，斜向45°照射時雙站RCS明顯小于正入射時雙站RCS，一定程度上反映了斜向入射情況下電磁泄漏要小于正入射.將360個雙站RCS值進行算術平均，RCS均值表如表1所示.

表1 雙站RCS全向均值統計 dBsm

由表1可知，入射角45°下RCS遠低于入射角0°，ft=1.5時 RCS遠低于 ft=10.可見，ft與入射角對電磁泄漏影響很大，且由表1可知其差別甚至接近30 dB.

3.2 時間網格因子的影響

下面研究ft對電磁泄漏的影響.計算過程中取入射角45°，ft取值從 1.5 到10，步長 1.5，其他參數同前.RCS均值隨ft的變化曲線如圖7所示.

圖7 RCS均值隨ft變化

由圖7可知，ft從1.5 變化到2.5，RCS均值迅速增加，而ft超過2.5之后，RCS均值的變化區域平緩，逐漸趨于一個固定的值，最后維持在約－18 dB左右.ft取值接近穩定性條件，可使遠場RCS降低約30 dB.

3.3 入射角度影響

二維情形下電磁泄漏與入射角度也密切相關.下面研究隨入射角度變化的影響，網格劃分同上.由于對稱性，入射角取為0°～45°就能夠代表所有的情形.入射角步長取為5°.

用遠場外推方法計算雙站RCS，雙站RCS均值隨入射角的變化如圖8所示.

圖8 雙站RCS均值隨入射角度的變化

由圖8可見，雙站RCS均值隨入射角度增大而減小.在入射角較小時，雙站RCS均值隨入射角改變比較緩慢，入射角大于15°后，雙站RCS均值改變逐漸加快.當入射角為45°時，達到最小值，約－37 dB.不同的入射角雙站RCS均值差別可達30 dB.

4 結論

本文用不同的衡量方法對一維、二維情形下連接邊界導致的電磁泄漏進行了研究，主要結論如下:

1)為減小電磁泄漏，時間步長在滿足穩定性要求的前提下，應越小越好.在二維情形下，時間步長過大引起的電磁泄漏可比時間步長接近穩定性條件時大30 dB.電磁泄漏隨時間步長變化的趨勢先慢后快，時間網格因子小于2時，基本呈線性關系.時間網格因子大于3后，電磁泄漏隨時間步長的增加緩慢變化，趨于穩定

2)入射角對電磁泄漏也有很大影響.當入射波斜入射到連接邊界時，不同的入射角對電磁泄漏也有很大影響.當入射角與計算邊界成斜角時，電磁泄漏較小.對于正方形區域，入射角度45°時電磁泄漏最小.垂直入射和斜入射時，電磁泄漏隨入射角度變化較為緩和，其他角度變化較劇烈.

References)

［1］鄭奎松，葛德彪，魏兵.導彈目標的FDTD建模與RCS計算［J］.系統工程與電子技術，2004，26(7):896－899 Zheng Kuisong，Ge Debiao，Wei Bing.FDTD modeling of missile target and RCS computation［J］.Systems Engineering and Elec－tronics，2004，26(7):896－899(in Chinese)

［2］李軍，武振波，武哲.穩定因子對FDTD數值計算的影響［J］.北京航空航天大學學報，2004，30(1):70－73 Li Jun，Wu Zhenbo，Wu Zhe.Effects of stability factor on FDTD numerical calculation［J］.Journal of Beijing University of Aeronautics and Astronautics，2004，30(1):70－73(in Chinese)

［3］黃志祥，沙威，吳先良，等.辛FDTD 算法［J］.系統工程與電子技術，2009，31(2):456－458 Huang Zhixiang，Sha Wei，Wu Xianliang，et al.Scheme of symplectic FDTD［J］.Systems Engineering and Electronics，2009，31(2):456－458(in Chinese)

［4］胡曉娟，盧兆林，葛德彪，等.基于三角面元的涂層目標FDTD共形網格生成技術［J］.系統工程與電子技術，2010，32(9):1884－1888 Hu Xiaojuan，Lu Zhaolin，Ge Debiao，et al.Conformal FDTD mesh－generating scheme for coated targets based on trianglepatch［J］.Systems Engineering and Electronics，2010，32(9):1884－1888(in Chinese)

［5］葛德彪，閆玉波.電磁波時域有限差分方法［M］.2版.西安:西安電子科技大學出版社，2005 Ge Debiao，Yan Yubo.Finite－difference time－domain method for electromagnetic waves［M］.2nd ed.Xi′an:Xidian University Press，2005(in Chinese)

［6］王長清.現代計算電磁學基礎［M］.北京:北京大學出版社，2004 Wang Changqing.Basic of modern computational electromagnetics［M］.Beijing:Peking University Press，2004(in Chinese)

［7］姜彥南，葛德彪，魏兵.時域有限差分并行算法中的吸收邊界研究［J］.系統工程與電子技術，2008，30(9):1636－1640 Jiang Yannan，Ge Debiao，Wei Bing.Study on absorbing boundary condition in parallel FDTD algorithm［J］.Systems Engineering and Electronics，2008，30(9):1636－1640(in Chinese)

［8］Wang Hui，Huang Zhixiang，Wu Xianliang，et al.Perfect plane wave injection into 3D FDTD(2，4)scheme［C］//Cross Strait Quad－Regional Radio Science and Wireless Technology Conference(CSQRWC)Volume 1.Harbin:IEEE，2011:40－43

［9］Hadi M F.A versatile split－field 1－D propagator for perfect FDTD plane wave injection［J］.IEEE Transactions on Antennas and Propagation，2009，57(9):2691－2697

［10］Huang Z，Pan G，Pan H K.Perfect plane wave injection for Crank－Nicholson time－domain method［J］.Microwaves，Antennas ＆ Propagation，IET，2010，4(11):1855－1862