姬金祖

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

黃大慶

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

黃沛霖 魯振毅

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

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

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

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

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

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

1 FDTD計算區(qū)域劃分

計算目標(biāo)電磁散射過程中，F(xiàn)DTD計算區(qū)域一般劃分為總場區(qū)和散射場區(qū)，在散射場區(qū)還設(shè)置輻射邊界，用于將近場外推至遠場，進一步計算RCS.FDTD計算區(qū)域劃分如圖1所示.

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

圖1 FDTD計算區(qū)域劃分

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

2 一維情形電磁泄漏

2.1 差分方程

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

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

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

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

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

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

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

2.2 時間網(wǎng)格因子影響

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

圖2 不同時間網(wǎng)格因子電場分布

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

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

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

圖3 電磁泄漏能量隨時間網(wǎng)格因子的變化

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

圖4 電磁泄漏隨時間網(wǎng)格因子的變化

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

可見，在一維FDTD中，時間網(wǎng)格因子在滿足穩(wěn)定性條件的基礎(chǔ)上，取值應(yīng)盡量小，即接近于1.

3 二維情形電磁泄漏

3.1 二維情形電磁泄漏比較

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

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

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

圖5 入射角度示意圖

圖6 無目標(biāo)下遠場雙站RCS

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

表1 雙站RCS全向均值統(tǒng)計 dBsm

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

3.2 時間網(wǎng)格因子的影響

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

圖7 RCS均值隨ft變化

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

3.3 入射角度影響

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

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

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

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

4 結(jié)論

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

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

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

References)

［1］鄭奎松，葛德彪，魏兵.導(dǎo)彈目標(biāo)的FDTD建模與RCS計算［J］.系統(tǒng)工程與電子技術(shù)，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］李軍，武振波，武哲.穩(wěn)定因子對FDTD數(shù)值計算的影響［J］.北京航空航天大學(xué)學(xué)報，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］.系統(tǒng)工程與電子技術(shù)，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］胡曉娟，盧兆林，葛德彪，等.基于三角面元的涂層目標(biāo)FDTD共形網(wǎng)格生成技術(shù)［J］.系統(tǒng)工程與電子技術(shù)，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版.西安:西安電子科技大學(xué)出版社，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］王長清.現(xiàn)代計算電磁學(xué)基礎(chǔ)［M］.北京:北京大學(xué)出版社，2004 Wang Changqing.Basic of modern computational electromagnetics［M］.Beijing:Peking University Press，2004(in Chinese)

［7］姜彥南，葛德彪，魏兵.時域有限差分并行算法中的吸收邊界研究［J］.系統(tǒng)工程與電子技術(shù)，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