許勇，黃勇，孫俊峰

(中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000)

0 引 言

復(fù)雜外形目標的電磁散射具有復(fù)雜散射機理，通常包括鏡面散射、邊緣繞射、爬行波、表面波、幾何非連續(xù)結(jié)構(gòu)散射等復(fù)雜電磁現(xiàn)象，準確模擬非常復(fù)雜和困難。另一方面，發(fā)展快速高效、精度高的軍用復(fù)雜外形飛行器的電磁散射計算方法為氣動/隱身優(yōu)化設(shè)計所迫切需要。傳統(tǒng)的高頻漸進方法和部件分解法等雖有方法通用、快捷，易于目標建模和能捕獲主要散射特征的優(yōu)點，但也有計算誤差較大的缺點。目前國內(nèi)外發(fā)展的電磁散射數(shù)值方法主要包括兩類：一類是求解電流積分方程，典型的例如多層快多極子(MLFMA)方法[1-4]；另一類是求解微分方程的FDTD方法[5-7]以及有限元(FEM)方法。積分方程方法采用格林函數(shù)避免電磁波在空間傳播的耗散、色散誤差，但會帶來稠密的耦合系數(shù)矩陣；微分方程方法直接計算電磁場，涉及場的時空計算、傳播和累積誤差，但也有適用問題廣，易于編程的優(yōu)點。

電磁學麥克斯韋方程組和流體力學無粘流歐拉方程都是有實特征值的雙曲型偏微分方程組，相同的數(shù)學特性允許采用同樣的偏微分數(shù)值算法。時域有限體積法(FVTD)[8-11]廣泛應(yīng)用于CFD工程中，該方法采用貼體曲線坐標系，避免了傳統(tǒng)FDTD[5]中笛卡爾網(wǎng)格帶來的階梯效應(yīng)誤差，不同于FDTD二階中心差分格式和利用電磁場量時/空錯置來提供人工粘性，F(xiàn)VTD利用迎風格式和一定網(wǎng)格密度來降低數(shù)值計算中的耗散和色散誤差，其計算機存儲與未知量數(shù)目同數(shù)量級，此時間推進方法能相容地模擬散射、多重散射、孔穿透、腔激勵等復(fù)雜現(xiàn)象而不需特殊處理。國外，在前期導(dǎo)電體電磁場計算的基礎(chǔ)上，F(xiàn)VTD目前主要發(fā)展在應(yīng)用層面，例如集成電路設(shè)計，F(xiàn)VTD和高頻混合方法計算電大目標電磁散射。D.K.Firsov等[12]研究了FVTD和積分方程(IE)結(jié)合節(jié)省計算空間和網(wǎng)格量；A.Chatterjee[13]研究了代數(shù)多重網(wǎng)格技術(shù)和FVTD方法結(jié)合來有效模擬線性電磁波傳播及保證高階空間離散精度。國內(nèi)，前期主要的FVTD研究工作側(cè)重于完全導(dǎo)電體電磁散射計算和多學科優(yōu)化應(yīng)用[14-16]。2014年，聶在平等[17]提出：超電大目標、介質(zhì)/導(dǎo)體混合結(jié)構(gòu)、腔體以及多尺度電磁散射問題面臨工程挑戰(zhàn)，目前基于CFD的電磁解算器仍有所欠缺。

為此，本文擴展FVTD方法在目標電磁散射問題中的應(yīng)用，首先，發(fā)展寬帶信號入射電磁波FVTD方法，在一個非定常過程計算多個頻率電磁散射；然后，發(fā)展特殊通量計算方法計算介質(zhì)/金屬目標電磁散射問題；最后，發(fā)展并行FVTD方法數(shù)值計算電大尺寸目標電磁散射問題。

1 數(shù)值方法

任何電磁問題的電磁場解都滿足如下時變麥克斯韋方程組：

(1)

(2)

(3)

(4)

式中：B為磁感應(yīng)強度矢量；H為磁場強度；D為電位移矢量；E為電場強度矢量；J為自由電流密度；ρ為自由電荷密度。

無源情況下，例如在自由空間中，J=0，ρ=0。

時變麥克斯韋方程組在直角坐標系下的守恒形式為

(5)

其中，

對于復(fù)雜外形物體，經(jīng)坐標變換：

得到曲線坐標系下的麥克斯韋方程組守恒形式為

(6)

其中，

(7)

式中：V為坐標變換的雅可比矩陣行列式值。

有限體積法的特點是能保持整個網(wǎng)格空間的通量守恒，對守恒方程(6)在每個網(wǎng)格單元作時間、空間積分，可以得到方程的數(shù)值離散化形式：

(8)

有限體積法的空間精度體現(xiàn)在能否精確模擬原變量Q在網(wǎng)格單元分界面處的狀態(tài)變量，以得到相應(yīng)精確的分界面流通量F。流通量的計算采用Steger-Warming分裂，也可通過求分界面處的黎曼(Riemann)解等方法得到。

(9)

式(9)中的單元邊界左右流通量可統(tǒng)一寫為

(10)

(11)

式中：k=ξ,η,ζ；S,S-為相似矩陣；Λ+,Λ-分別為正負特征值構(gòu)成的對角矩陣；QL,QR分別為分界面處左右狀態(tài)變量，可采用MUSCL格式得到最高三階精度。

(12)

(13)

完全導(dǎo)電壁面反射邊界條件，由電磁理論為

(14)

(15)

式中：E,B為總場。

式(14)～式(15)是不完備的，未提供電場垂直表面、磁場切向于表面信息，補充近似邊界條件：

(16)

(17)

本文借鑒CFD中固壁邊界處理方式，發(fā)展了三種可行的完全導(dǎo)電體(PEC)邊界條件，分別是1、2階外插邊界條件以及虛擬像點方式，虛擬像點方式有通用性，有利于數(shù)據(jù)交換和并行通信。

截斷網(wǎng)格外邊界采用輻射邊界條件或吸收邊界條件，以降低邊界反射回波帶來的畸變。本文采用相容條件：

(18)

時間計算方面采用與常微分方程相似的龍格-庫塔方法：

Qn+k/m=Qn-λαkR(Qn+(k-1)/m) (k=1,m)

(19)

其中，

式中：m為龍格-庫塔法的步數(shù)，本文中m=4；R為方程的殘差。

2 電磁散射特性計算

2.1 寬帶脈沖電磁波入射

時域方法的一個重要優(yōu)點是在寬帶脈沖電磁波入射情況下，利用非定常計算和傅里葉變換，在一個計算狀態(tài)中可以獲得多個頻率目標電磁散射特性。

入射信號采用高斯脈沖：

(20)

頻譜強度由對應(yīng)的非周期信號傅里葉變換計算：

(21)

本文實際計算中從5倍半周期Tm信號，即5×Tm開始，整個網(wǎng)格密度由待求最高頻(最短波長)決定，如此低頻信號網(wǎng)格密度自然滿足精度要求，停止計算條件為積分面上最大電磁場幅度小于-18 dB。

以高斯型寬帶脈沖電磁波入射情況下的金屬球雷達散射截面計算為例，網(wǎng)格為49×121×61。輸入電磁信號如圖1所示，頻譜分布如圖2所示，金屬球?qū)拵盘栯姶派⑸渲须姶艌鰰r間歷程和3個頻率對應(yīng)雙站RCS分布如圖3所示，可以看出：與對應(yīng)解析解吻合很好。

圖1 入射高斯型脈沖電磁波信號Fig.1 Incident electromagnetic wave of Gaussian pulse

圖2 頻譜分布Fig.2 Spectrum of frequency domain

(a) 脈沖散射場時間歷程

(b) ka=1

(c) ka=5

(d) ka=10圖3 寬帶脈沖波入射金屬球雷達截面(RCS)計算Fig.3 Bistatic RCS profiles of incident wideband signal with different frequency

2.2 介質(zhì)/金屬目標電磁散射計算

介質(zhì)通常指不同于自由空間，其普遍形式是具有復(fù)數(shù)型介電常數(shù)和磁化率，實部引起電磁波折射，虛部引發(fā)電磁波耗散。采用散射電磁場守恒形式，電導(dǎo)率和磁導(dǎo)率分別為：σe=ωεi，σm=ωμi。利用介質(zhì)電磁參數(shù)間斷處電位移矢量、磁感應(yīng)強度切向分量連續(xù)插值和構(gòu)建通量。

(22)

(23)

一類不同磁化率介質(zhì)覆蓋金屬球的雙站RCS計算結(jié)果如圖4所示。計算條件：ka=6.28，介質(zhì)厚度：d=λd/30 。

(a) E平面

從圖4可以看出：介質(zhì)虛部引發(fā)電磁能量耗散，相應(yīng)降低散射電磁波能量，從而帶來RCS降低，這也是等離子體隱身機理的驗證。

2.3 電大尺寸金屬目標電磁散射計算

FVTD方法直接求解麥克斯韋方程組，是全波數(shù)值方法，適用于從低頻到高頻全范圍，但由于空間網(wǎng)格數(shù)量與頻率平方成正比，對每個波長如取大于15個網(wǎng)格點情況下，戰(zhàn)斗機X波段網(wǎng)格數(shù)量達到數(shù)十億規(guī)模。對超電大尺寸目標必須輔以并行算法，因涉及電磁波在三維計算空間傳播，其完全導(dǎo)電體電磁散射計算效率比不上多層快多極子方法(MLFMA)。

高頻電大尺寸目標電磁散射問題中，三維空間網(wǎng)格的大幅增長使得計算量十分龐大，并行計算勢在必行，因此構(gòu)建多進程并行平臺，包括：①網(wǎng)格多進程分割和負載平衡，②程序并行化處理，采用MPI接口進行通信。文獻[11]詳細介紹了該并行解算器的驗證和相關(guān)細節(jié)。本文驗證計算算例為美國EMCC標模。歸一化表面誘導(dǎo)電流等值線云圖如圖5所示，顯示電磁反射強度區(qū)域分布，截斷錐臺單站RCS計算與測量比較如圖6所示。高度200 mm,底部直徑200、100 mm,計算頻率7 GHz,平行極化，計算驗證FVTD模擬邊緣繞射和爬行波的能力。使用多塊結(jié)構(gòu)網(wǎng)格，包含257萬網(wǎng)格點，采用100個進程，由于并行FVTD程序結(jié)構(gòu)同于CFD解算器，其并行效率亦同于CFD中并行有限體積法流場解算器。

圖5 歸一化表面誘導(dǎo)電流分布Fig.5 Contuor of normalized induced surface current

圖6 截斷錐臺單站RCS計算與測量比較(f=7 GHz)Fig.6 Backscattering RCS of truncated cone compared with measurement(f=7 GHz)

3 結(jié) 論

(1) 針對寬帶脈沖入射電磁波，研究通過非周期傅里葉變換和非定常FVTD計算，在一個計算狀態(tài)中獲得多個頻率電磁散射特性，金屬球雙站RCS與解析解驗證比較，吻合良好。

(2) 對于涂覆吸波材料類的介質(zhì)/導(dǎo)體電磁散射FVTD模擬，導(dǎo)電率在分界面存在間斷，相應(yīng)電磁場在介質(zhì)參數(shù)突變也存在間斷，利用物理邊界條件，研究間斷點明確的通量和特殊的插值方式，計算算例中與Mie級數(shù)解誤差小于1 dB。

(3) 針對高頻、電大尺寸目標研究了相應(yīng)并行算法，包括網(wǎng)格多進程分割和負載平衡、程序并行化處理及MPI接口通信，成功計算高頻目標電磁散射，并與暗室測量結(jié)果吻合良好，誤差低于1.5 dBsm。