CFD技術在目標電磁特性計算中的應用

許勇，黃勇，孫俊峰

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

0 引 言

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

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

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

1 數值方法

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

(1)

(2)

(3)

(4)

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

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

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

(5)

其中，

對于復雜外形物體，經坐標變換：

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

(6)

其中，

(7)

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

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

(8)

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

(9)

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

(10)

(11)

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

(12)

(13)

完全導電壁面反射邊界條件，由電磁理論為

(14)

(15)

式中：E,B為總場。

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

(16)

(17)

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

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

(18)

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

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

(19)

其中，

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

2 電磁散射特性計算

2.1 寬帶脈沖電磁波入射

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

入射信號采用高斯脈沖：

(20)

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

(21)

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

以高斯型寬帶脈沖電磁波入射情況下的金屬球雷達散射截面計算為例，網格為49×121×61。輸入電磁信號如圖1所示，頻譜分布如圖2所示，金屬球寬帶信號電磁散射中電磁場時間歷程和3個頻率對應雙站RCS分布如圖3所示，可以看出：與對應解析解吻合很好。

圖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 介質/金屬目標電磁散射計算

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

(22)

(23)

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

(a) E平面

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

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

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

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

圖5 歸一化表面誘導電流分布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 結 論

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

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

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