基于CFD無(wú)網(wǎng)格算法的飛行器電磁隱身特性模擬

高煜堃，陳紅全，王 彪，吳 浩，胡曉磊，宋強(qiáng)強(qiáng)

(1.安徽工業(yè)大學(xué)機(jī)械工程學(xué)院，安徽馬鞍山243002；2.南京航空航天大學(xué)航空宇航學(xué)院，江蘇南京 210016)

關(guān)于現(xiàn)代先進(jìn)軍用飛行器電磁隱身特性的研究一直是學(xué)術(shù)界和工程界關(guān)注的熱點(diǎn)問(wèn)題。雷達(dá)散射截面(radar cross section，RCS)作為評(píng)價(jià)飛行器電磁隱身特性的一項(xiàng)重要指標(biāo)，通常通過(guò)求解麥克斯韋方程獲得。隨著計(jì)算機(jī)水平的不斷提高，越來(lái)越多的研究者開(kāi)始嘗試采用數(shù)值方法來(lái)求解麥克斯韋方程，獲取目標(biāo)RCS。傳統(tǒng)的數(shù)值方法有時(shí)域有限差分法[1]、間斷有限元法[2]、時(shí)域有限體積法[3-4]和時(shí)域有限元法[5]等，這些方法均依賴(lài)于網(wǎng)格單元離散計(jì)算域，受到網(wǎng)格化的約束。無(wú)網(wǎng)格方法打破了傳統(tǒng)網(wǎng)格方法受到的網(wǎng)格化約束，其計(jì)算域的離散只涉及布點(diǎn)，不需生成網(wǎng)格單元，實(shí)施靈活，適合于處理多體干擾及多部件干擾等復(fù)雜情形問(wèn)題。

目前，已有研究者嘗試采用無(wú)網(wǎng)格方法進(jìn)行電磁場(chǎng)數(shù)值計(jì)算，其中典型的有無(wú)單元Galerkin方法[6]和徑向基函數(shù)的無(wú)網(wǎng)格方法[7]。上述兩類(lèi)無(wú)網(wǎng)格求解方法均與基函數(shù)(或形函數(shù))的選取相關(guān)，基函數(shù)的選取不僅會(huì)影響矩陣求逆運(yùn)算中涉及的逆矩陣質(zhì)量，進(jìn)而影響算法的收斂性；而且會(huì)影響到邊界條件的處理，有時(shí)邊界條件需采用特殊的方法強(qiáng)制給定。在計(jì)算流體力學(xué)(computational fluid dynamics，CFD)領(lǐng)域，近年來(lái)出現(xiàn)了一種用于求解歐拉方程的無(wú)網(wǎng)格方法[8],其空間離散不涉及基函數(shù)的選取，避免了逆矩陣對(duì)算法收斂性的影響，邊界條件的施加也相對(duì)簡(jiǎn)單，通用性較好?；诖?，筆者借鑒CFD無(wú)網(wǎng)格方法的思想，提出一種基于CFD的時(shí)域無(wú)網(wǎng)格算法，用于求解麥克斯韋方程，采用Steger-Warming通量分裂進(jìn)行通量計(jì)算，且采用本文算法對(duì)文獻(xiàn)[9]中的平板飛機(jī)模型進(jìn)行電磁隱身特性模擬。

1 基于CFD的時(shí)域無(wú)網(wǎng)格算法

1.1 控制方程

對(duì)于橫磁(transverse magnetic，TM)波，在直角坐標(biāo)系中，守恒型量綱為一的時(shí)域麥克斯韋方程[10]可寫(xiě)為

其中：W=[εEzμHxμHy]T；F1=[-Hy0 -Ez]T；F2=[HxEz0]T；S=[-σEz-σmHx-σmHy]T；Ez為電場(chǎng)強(qiáng)度E在z方向上的分量；Hx和Hy分別為磁場(chǎng)強(qiáng)度H在x和y方向上的分量；ε為介電常數(shù)；μ為磁導(dǎo)率；σ為電導(dǎo)率；σm為磁阻率。

1.2 計(jì)算區(qū)域布點(diǎn)及點(diǎn)云生成

與傳統(tǒng)的網(wǎng)格方法不同，無(wú)網(wǎng)格方法計(jì)算區(qū)域的離散只涉及布點(diǎn)而不需生成網(wǎng)格單元，通常情況下，既可直接采用網(wǎng)格點(diǎn)對(duì)計(jì)算域進(jìn)行離散，也可根據(jù)需要進(jìn)行布點(diǎn)離散計(jì)算域。對(duì)計(jì)算域布點(diǎn)離散后需生成局部點(diǎn)云結(jié)構(gòu)。以點(diǎn)云Ci為例(見(jiàn)圖1)，文獻(xiàn)[11]中給出了點(diǎn)云Ci的具體生成方法。其中點(diǎn)i為點(diǎn)云Ci的中心點(diǎn)，點(diǎn)1～6為點(diǎn)云Ci的衛(wèi)星點(diǎn)。

圖1 點(diǎn)云Ci示意圖Fig.1 Schematic diagram of cloud of pointsCi

1.3 空間導(dǎo)數(shù)逼近

在點(diǎn)云Ci上，假設(shè)中心點(diǎn)i附近的函數(shù)分布f=f(x,y)和連續(xù)可微，那么f可用i處的函數(shù)值fi=f(xi,yi)通過(guò)泰勒級(jí)數(shù)展開(kāi)逼近

其中：h=x-xi；h=y-yi；ai(i=1,2)為函數(shù)f在中心點(diǎn)i處的空間導(dǎo)數(shù)。對(duì)于線性逼近，式(2)函數(shù)的近似值可寫(xiě)為fˉ=fi+a1h+a2l。對(duì)于時(shí)域無(wú)網(wǎng)格算法，衛(wèi)星點(diǎn)k(k=1,…,M)處的函數(shù)值為已知，記作fk。為使衛(wèi)星點(diǎn)處函數(shù)值的線性逼近總體誤差取到極小值，則空間導(dǎo)數(shù)滿(mǎn)足下式[11]

那么線性逼近函數(shù)可整理為f(x,y)=fi+∑αk(fk-fi)h+βk(fk-fi)l，其中系數(shù)αk和βk僅與離散點(diǎn)坐標(biāo)相關(guān)，在時(shí)間推進(jìn)計(jì)算前可一次求出。因此，函數(shù)f在中心點(diǎn)i處的空間導(dǎo)數(shù)可逼近為

空間導(dǎo)數(shù)也可用中心點(diǎn)與衛(wèi)星點(diǎn)連線中點(diǎn)處的函數(shù)值進(jìn)行逼近

其中系數(shù)αik和βik在時(shí)間推進(jìn)計(jì)算前也可一次求出。

1.4 通量運(yùn)算

在點(diǎn)云Ci上，運(yùn)用式(5)，則中心點(diǎn)i處的通量項(xiàng)可寫(xiě)為

由于求和項(xiàng) ∑(αikF1k+βikF2i)為已知(系數(shù)αik和βik在時(shí)間推進(jìn)計(jì)算前已經(jīng)計(jì)算得到，F(xiàn)1i和F2i為計(jì)算點(diǎn)處當(dāng)前的已知值)，在中心點(diǎn)與每一個(gè)衛(wèi)星點(diǎn)連線的中點(diǎn)處建立一個(gè)虛擬界面，如圖2。并定義一個(gè)數(shù)值通量Qik=ξikF1(Wik)+ηikF2(Wik)，其中定義矢量dik=(αik，βik）,那么式(6)可寫(xiě)為

圖2 點(diǎn)云Ci中心點(diǎn)與衛(wèi)星點(diǎn)連線中點(diǎn)處引入的虛擬界面Fig.2 Virtual interface between the central and each satellite point on the point cloud ofCi

為了計(jì)算Qik，借鑒CFD的做法，采用Steger-Warming通量分裂方法[4]分裂Qik的雅克比矩陣Tik=?Qik/?Wik，那么Qik可寫(xiě)成分裂形式其中為T(mén)ik的分裂矩陣，可由線性函數(shù)重構(gòu)得到U+=Ui+0.5?Ui?rik，U-=Uk+0.5?Uk?rik，其中rik=(xk-xi，yk-yi)，這里U表示W(wǎng)的任意分量，?Ui和?Uk可以由式(4)計(jì)算得到。

1.5 時(shí)間推進(jìn)與邊界條件

空間離散后，在點(diǎn)云Ci上，麥克斯韋方程的半離散形式可寫(xiě)為

其中Ri為計(jì)算點(diǎn)i處的殘差。然后按照四步Runge-Kutta方法[4]對(duì)式(8)進(jìn)行時(shí)間推進(jìn)求解。文中物面采用良導(dǎo)體邊界條件，截?cái)噙吔绮捎猛耆ヅ鋵舆吔鐥l件，具體參數(shù)取值同文獻(xiàn)[4]。

圖3 入射波示意圖Fig.3 Schematic diagram of incident wave

2 算例及結(jié)果分析

算法采用fortran編程實(shí)現(xiàn)，并采用該程序?qū)ξ墨I(xiàn)[9]中飛機(jī)模型的隱身特性進(jìn)行分析。圖3為沿著k′方向傳播的TM平面波示意圖，定義k′軸與x軸之間的夾角為入射角φ，歸一化后的入射波分量可表示為

其中k′=xcosφ×ysinφ。算例涉及的計(jì)算域中x和y坐標(biāo)均為以入射波波長(zhǎng)λ為特征長(zhǎng)度、量綱為一的空間位置。

2.1 二維圓柱散射數(shù)值模擬

圖4 二維圓柱的散射場(chǎng)分布Fig.4 Distribution of the scattered fields of 2D cylinders

選用二維圓柱算例對(duì)本文算法進(jìn)行驗(yàn)證。數(shù)值模擬時(shí)，取圓柱半徑r=0.5λ，置于10λ×10λ的計(jì)算域中，總布點(diǎn)數(shù)為14 458，采用沿x軸方向的TM波入射(對(duì)應(yīng)φ=0°)。計(jì)算得到的圓柱散射場(chǎng)分布和雙站RCS分布分別見(jiàn)圖4，5。由圖4可知，當(dāng)電磁波沿x軸方向照射時(shí)，雙站角0°方向的散射場(chǎng)最強(qiáng)，這與計(jì)算得到的雙站RCS分布(見(jiàn)圖5)一致，RCS最大值出現(xiàn)在雙站角0°方向。由圖5可知，基于本文算法計(jì)算得到的雙站RCS分布在整個(gè)雙站角內(nèi)都能與級(jí)數(shù)解[12]吻合。圖6為二維圓柱的殘值收斂歷程。由圖6可知，當(dāng)?shù)芷跒?6時(shí)，本文算法的迭代平均殘值下降到1.0×107以下，收斂速度良好。

圖5 二維圓柱的雙站RCS分布Fig.5 Distribution of the bistatic RCS of 2D cylinders

圖6 二維圓柱的殘值收斂歷程Fig.6 Convergence process of the value of 2D cylinders

2.2 平板飛機(jī)模型隱身特性模擬

為驗(yàn)證本文算法處理多體及多部件干擾問(wèn)題的能力，選用文獻(xiàn)[9]中的平板飛機(jī)模型進(jìn)行電磁隱身特性模擬。數(shù)值模擬時(shí)，將兩架相同的、長(zhǎng)度均為6λ的平板飛機(jī)模型并排放置在大小為24λ×24λ的計(jì)算域中，見(jiàn)圖7，總布點(diǎn)數(shù)為46 917。為研究電磁波從不同方向照射時(shí)飛機(jī)模型的隱身特性，取TM波入射角φ=0°～90°范圍內(nèi)(對(duì)應(yīng)電磁波從飛機(jī)模型的側(cè)前方照射)每間隔5°的19種情形，從RCS平均值、最大值以及大于門(mén)檻值的概率3個(gè)指標(biāo)來(lái)分析飛機(jī)模型的隱身特性，根據(jù)文獻(xiàn)[9]設(shè)置RCS門(mén)檻值分析RCS隱身特性，文中將該值設(shè)置為0 dB，對(duì)應(yīng)計(jì)算得到的雙站RCS結(jié)果見(jiàn)表1。一般來(lái)說(shuō)，RCS平均值越小越好，RCS最大值越小越好，RCS大于門(mén)檻值的概率越低越好。

圖7 飛機(jī)模型計(jì)算域布點(diǎn)離散示意圖Fig.7 Points distribution in the computational domain for the aircraft models

由表1第四列可知，RCS最大值出現(xiàn)的角度與入射角相同，表明在入射角方向上的散射最強(qiáng)。由表1第二列可知：φ=0°時(shí)，RCS平均值最大，為7.427 dB；φ=5°，15°時(shí)，RCS平均值較大；φ=50°時(shí)，RCS平均值最小，為4.468 dB；φ=65°，90°時(shí)，RCS平均值較小。由表1第三列可知：φ=45°時(shí)，RCS最大值最大為28.116 dB；φ=30°，35°時(shí)，RCS最大值較大；φ=90°時(shí)，RCS最大值最小，為26.006 dB；φ=80°，85°時(shí)，RCS最大值較小。由表1第五列可知：φ=0°，5°時(shí)，RCS大于門(mén)檻值的概率最大，為88.1%；φ=50°時(shí)，RCS大于門(mén)檻值的概率最小，為73.3%；φ=65°，90°時(shí)，RCS大于門(mén)檻值的概率較小。綜合來(lái)看：φ=0°時(shí)(對(duì)應(yīng)電磁波從正前方照射)，RCS平均值和大于門(mén)檻值的概率均最大，隱身特性較差；φ=90°時(shí)(對(duì)應(yīng)電磁波從側(cè)向照射)，RCS平均值、最大值以及大于門(mén)檻值的概率均較小，隱身特性較好。

表1 飛機(jī)模型不同入射角情形對(duì)應(yīng)的雙站RCS分析Tab.1 Analysis of the bistatic RCS of aircraft models with different incident angles

圖8 飛機(jī)模型散射場(chǎng)分布Fig.8 Distribution of the scattered fields of aircraft models

圖8，9分別為飛機(jī)模型在φ=0°和φ=90°情形下散射場(chǎng)分布和雙站RCS分布。由圖8可知：當(dāng)電磁波從正前方照射飛機(jī)模型時(shí)(φ=0°)，尾翼被機(jī)翼遮擋，機(jī)翼之前的外形布局部件間二面角均大于90°,可有效回避散射波的疊加，表現(xiàn)為各部件產(chǎn)生的散射波在散射場(chǎng)中有序發(fā)散傳播(圖8(a))，最強(qiáng)散射出現(xiàn)在雙站角0°方向，這與圖9中的RCS分布一致；當(dāng)電磁波從側(cè)向照射飛機(jī)模型時(shí)(φ=90°)，機(jī)頭與機(jī)身、機(jī)身與機(jī)翼以及機(jī)身與尾翼之間都存在一定的電磁散射干擾，在雙站角180°～360°范圍內(nèi)最強(qiáng)散射方向上的散射強(qiáng)度被削弱，故在整個(gè)雙站角范圍內(nèi)最強(qiáng)散射只出現(xiàn)在90°方向上(圖8(b))，這也與圖9中的RCS分布吻合。

圖9 飛機(jī)模型雙站RCS分布Fig.9 Distribution of the bistatic RCS of aircraft models

3 結(jié) 論

借鑒CFD中Steger-Warming通量分裂方法，提出基于CFD的時(shí)域無(wú)網(wǎng)格算法，并用于分析飛行器的電磁隱身特性，結(jié)論如下：

1)采用本文算法計(jì)算得到的二維圓柱雙站RCS能與級(jí)數(shù)解吻合；

2)從RCS平均值、最大值以及大于門(mén)檻值的概率三方面綜合分析飛機(jī)模型的隱身特性，當(dāng)入射角φ=0°時(shí)(對(duì)應(yīng)電磁波從正前方照射)飛機(jī)隱身特性較差，而當(dāng)入射角φ=90°時(shí)(對(duì)應(yīng)電磁波從側(cè)向照射)飛機(jī)隱身特性較好，這與散射場(chǎng)的疊加作用以及飛機(jī)外形等因素有關(guān)；

3)本文算法基于無(wú)網(wǎng)格點(diǎn)云構(gòu)造，適合于處理多體及多部件目標(biāo)等電磁散射干擾問(wèn)題。