基于CFD無網格算法的飛行器電磁隱身特性模擬

高煜堃，陳紅全，王 彪，吳 浩，胡曉磊，宋強強

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

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

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

1 基于CFD的時域無網格算法

1.1 控制方程

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

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

1.2 計算區域布點及點云生成

與傳統的網格方法不同，無網格方法計算區域的離散只涉及布點而不需生成網格單元，通常情況下，既可直接采用網格點對計算域進行離散，也可根據需要進行布點離散計算域。對計算域布點離散后需生成局部點云結構。以點云Ci為例(見圖1)，文獻[11]中給出了點云Ci的具體生成方法。其中點i為點云Ci的中心點，點1～6為點云Ci的衛星點。

圖1 點云Ci示意圖Fig.1 Schematic diagram of cloud of pointsCi

1.3 空間導數逼近

在點云Ci上，假設中心點i附近的函數分布f=f(x,y)和連續可微，那么f可用i處的函數值fi=f(xi,yi)通過泰勒級數展開逼近

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

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

空間導數也可用中心點與衛星點連線中點處的函數值進行逼近

其中系數αik和βik在時間推進計算前也可一次求出。

1.4 通量運算

在點云Ci上，運用式(5)，則中心點i處的通量項可寫為

由于求和項 ∑(αikF1k+βikF2i)為已知(系數αik和βik在時間推進計算前已經計算得到，F1i和F2i為計算點處當前的已知值)，在中心點與每一個衛星點連線的中點處建立一個虛擬界面，如圖2。并定義一個數值通量Qik=ξikF1(Wik)+ηikF2(Wik)，其中定義矢量dik=(αik，βik）,那么式(6)可寫為

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

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

1.5 時間推進與邊界條件

空間離散后，在點云Ci上，麥克斯韋方程的半離散形式可寫為

其中Ri為計算點i處的殘差。然后按照四步Runge-Kutta方法[4]對式(8)進行時間推進求解。文中物面采用良導體邊界條件，截斷邊界采用完全匹配層邊界條件，具體參數取值同文獻[4]。

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

2 算例及結果分析

算法采用fortran編程實現，并采用該程序對文獻[9]中飛機模型的隱身特性進行分析。圖3為沿著k′方向傳播的TM平面波示意圖，定義k′軸與x軸之間的夾角為入射角φ，歸一化后的入射波分量可表示為

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

2.1 二維圓柱散射數值模擬

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

選用二維圓柱算例對本文算法進行驗證。數值模擬時，取圓柱半徑r=0.5λ，置于10λ×10λ的計算域中，總布點數為14 458，采用沿x軸方向的TM波入射(對應φ=0°)。計算得到的圓柱散射場分布和雙站RCS分布分別見圖4，5。由圖4可知，當電磁波沿x軸方向照射時，雙站角0°方向的散射場最強，這與計算得到的雙站RCS分布(見圖5)一致，RCS最大值出現在雙站角0°方向。由圖5可知，基于本文算法計算得到的雙站RCS分布在整個雙站角內都能與級數解[12]吻合。圖6為二維圓柱的殘值收斂歷程。由圖6可知，當迭代周期為36時，本文算法的迭代平均殘值下降到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 平板飛機模型隱身特性模擬

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

圖7 飛機模型計算域布點離散示意圖Fig.7 Points distribution in the computational domain for the aircraft models

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

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

圖8 飛機模型散射場分布Fig.8 Distribution of the scattered fields of aircraft models

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

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

3 結 論

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

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

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

3)本文算法基于無網格點云構造，適合于處理多體及多部件目標等電磁散射干擾問題。