CFD技術(shù)在目標(biāo)電磁特性計(jì)算中的應(yīng)用

許勇，黃勇，孫俊峰

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000)

0 引 言

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

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

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

1 數(shù)值方法

任何電磁問(wèn)題的電磁場(chǎng)解都滿(mǎn)足如下時(shí)變麥克斯韋方程組：

(1)

(2)

(3)

(4)

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

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

時(shí)變麥克斯韋方程組在直角坐標(biāo)系下的守恒形式為

(5)

其中，

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

得到曲線(xiàn)坐標(biāo)系下的麥克斯韋方程組守恒形式為

(6)

其中，

(7)

式中：V為坐標(biāo)變換的雅可比矩陣行列式值。

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

(8)

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

(9)

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

(10)

(11)

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

(12)

(13)

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

(14)

(15)

式中：E,B為總場(chǎng)。

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

(16)

(17)

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

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

(18)

時(shí)間計(jì)算方面采用與常微分方程相似的龍格-庫(kù)塔方法：

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

(19)

其中，

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

2 電磁散射特性計(jì)算

2.1 寬帶脈沖電磁波入射

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

入射信號(hào)采用高斯脈沖：

(20)

頻譜強(qiáng)度由對(duì)應(yīng)的非周期信號(hào)傅里葉變換計(jì)算：

(21)

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

以高斯型寬帶脈沖電磁波入射情況下的金屬球雷達(dá)散射截面計(jì)算為例，網(wǎng)格為49×121×61。輸入電磁信號(hào)如圖1所示，頻譜分布如圖2所示，金屬球?qū)拵盘?hào)電磁散射中電磁場(chǎng)時(shí)間歷程和3個(gè)頻率對(duì)應(yīng)雙站RCS分布如圖3所示，可以看出：與對(duì)應(yīng)解析解吻合很好。

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

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

(a) 脈沖散射場(chǎng)時(shí)間歷程

(b) ka=1

(c) ka=5

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

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

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

(22)

(23)

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

(a) E平面

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

2.3 電大尺寸金屬目標(biāo)電磁散射計(jì)算

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

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

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

圖6 截?cái)噱F臺(tái)單站RCS計(jì)算與測(cè)量比較(f=7 GHz)Fig.6 Backscattering RCS of truncated cone compared with measurement(f=7 GHz)

3 結(jié) 論

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

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

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