等離子體鞘套電磁散射特性的場離散蒙特卡羅仿真

賈潔姝， 梁子長， 何鴻飛

（電磁散射重點實驗室，上海200438）

0 引言

飛行器超高速飛行時，與周圍大氣相互作用產(chǎn)生等離子體鞘套。鞘套中的帶電粒子發(fā)生碰撞會造成能量的損耗，對飛行器的電磁散射特性產(chǎn)生嚴重影響［1-3］。隨著飛行高度及馬赫數(shù)的變化，目標周圍所形成的繞流流場參數(shù)及電參數(shù)分布也發(fā)生改變，所呈現(xiàn)的目標雷達特征信息是不同的［4-5］。等離子體鞘套是影響探測超高速飛行器的重要因素。

超高速目標流場網(wǎng)格一般幾倍于目標尺寸，其網(wǎng)格數(shù)量較多。對等離子體中局部分層界面電磁波的計算，采用射線跟蹤方法需要進行射線分裂。對復(fù)雜的流場網(wǎng)格，該過程將十分復(fù)雜，計算效率也將下降。為解決這一問題，本文引入了蒙特卡羅方法（Monte Carlomethod），計算等離子體包覆目標電磁特性。

1 仿真建模

蒙特卡羅方法是一種以概率統(tǒng)計理論為指導(dǎo)的數(shù)值計算方法，在粒子輸運計算、熱力學(xué)、光學(xué)以及空氣動力學(xué)等領(lǐng)域應(yīng)用較為廣泛。由于蒙特卡羅方法是一種非確定性的算法，模擬對象需粒子化，在電磁散射計算領(lǐng)域，主要用于隨機環(huán)境的散射模擬，且主要基于能量離散電磁波粒子開展計算。

首先，以等離子體流場網(wǎng)格外表面為模擬電磁波粒子的入射面，在該入射面內(nèi)，電磁波粒子的起始位置呈隨機分布，粒子起始運動方向由電磁波源指向其起始位置。對每一電磁波粒子，假定在等離子體中的運動步長為Δl，則每前進一步，電磁波粒子位置變化Δζ及運動方向的變化Δθ分別如式（1）與式（2）所示。

式中：θ1為局部電磁波粒子前進方向υ與等離子體等效折射率變化梯度方向ζ的夾角；Δn為粒子向前行進一步對應(yīng)的等效折射率變化；n1為該位置處等離子體的等效折射率，在等離子體鞘套外表面處n1≈1。

然后，計算每前進一步電磁波粒子的吸收概率pa，如式（4）所示。

式中：a1為該位置處等離子體的等效吸收系數(shù)。利用計算機生成0至1間的偽隨機數(shù)，當該隨機數(shù)小于pa，則電磁波粒子被吸收。

再計算運動一步后電磁波粒子位置ζ1以及等效折射率n2，分別如式（5）與式（6）所示。

若n2＞0，電磁波粒子運動方向由式（7）計算。

式中：t為單位矢量，其方向為υ與ζ所在平面內(nèi)垂直于ζ的方向。

若n2≤0，則電磁波粒子運動方向發(fā)生全反射，如式（8）所示。

重復(fù)上述過程，直至電磁波粒子從等離子體網(wǎng)格外表面離開或被吸收。

重復(fù)以上過程，直至被吸收和被反射的電磁波粒子數(shù)與總?cè)肷淞Ｗ訑?shù)的比值趨于穩(wěn)定，一般要求該比值偏差小于10%。

等離子體包覆目標的RCS由被反射的電磁波粒子與入射粒子總數(shù)的比值計算，如式（9）所示。

式中：N為入射粒子總數(shù)，N越大，單個電磁波粒子表示的電場幅度越小，電場幅度在粒子傳播過程中保持不變；ΔEsi（θ）表示某觀測方向下接收的第i個散射粒子電場，其中包含路程相位和目標反射相位；M表示該方向散射粒子的總數(shù)。

與其它計算方法相比，該方法具有以下優(yōu)點：

a）適于復(fù)雜三維非均勻連續(xù)分布的等離子體及其包覆目標的電磁特性計算；

b）適于并行化計算，計算耗時與電磁波頻率相關(guān)性較小，可計算電大尺寸等離子體包覆目標的電磁特性；

c）可共用流場仿真網(wǎng)格，減少了計算網(wǎng)格轉(zhuǎn)換對電磁計算精度的影響。

2 收斂性分析

由于統(tǒng)計仿真方法存在漲落噪聲，MC方法計算等離子體目標的RCS需進行收斂性分析，即分析不同模擬粒子數(shù)對RCS結(jié)果的影響。這里以飛行馬赫數(shù)為15的鈍錐體5 GHz的RCS數(shù)據(jù)為例，按RCS峰值區(qū)域和谷值區(qū)域分別進行收斂性分析，如圖1和圖2所示。

圖1 鈍錐體RCS峰值區(qū)域的收斂性分析

在RCS峰值區(qū)域，當模擬粒子數(shù)由3萬增加至60萬時，等離子體覆蓋目標的RCS計算結(jié)果趨于穩(wěn)定，如圖1（b）所示。在RCS谷值區(qū)域，當模擬粒子數(shù)由1萬增加至60萬時，漲落噪聲略有降低，RCS計算結(jié)果趨于穩(wěn)定，如圖2（b）所示?？梢钥闯?，MC方法收斂性在RCS谷值和峰值區(qū)域有較大區(qū)別，在RCS谷值區(qū)域，粒子數(shù)的收斂性較慢，而在RCS峰值區(qū)域粒子數(shù)的收斂性十分快。

圖2 鈍錐體RCS谷值區(qū)域的收斂性分析

3 典型超高速目標模型RCS計算分析

采用上述的MC方法，對鈍錐體典型超高速目標計算其不同頻率、觀測角下的RCS，并對其中主要目標RCS衰減、增強等現(xiàn)象進行了分析。由于對稱性，流場網(wǎng)格僅提供一半網(wǎng)格，以下計算的RCS結(jié)果也僅考慮上半部分網(wǎng)格，實際目標的RCS應(yīng)增加約6 d B。

鈍錐體模型長為0.5 m，底面直徑0.33 m，如圖3所示。等離子體流場區(qū)域尺寸約為2.1 m×1.6 m×0.8 m。目標飛行馬赫數(shù)選取15和20，飛行高度為75 k m。以下分析不同頻率及方位角對RCS的影響。

主要計算參數(shù)：俯仰角為90°，即觀測面處于XOY面（即對稱面）內(nèi)；模擬入射粒子數(shù)為10萬；模擬粒子運動的最大步數(shù)為300。計算采用了6核并行，平均每頻點耗時約1.3 h。

圖3 鈍錐體模型

飛行馬赫數(shù)為15時，對比有無等離子體時鈍錐體RCS的變化，如圖4所示。頻率1 GHz在方位角20°附近，2 GHz在方位角5°附近，后向RCS出現(xiàn)無值情況，對應(yīng)完全的折射隱身現(xiàn)象。

同時，在頻率2 GHz、方位角20°以及頻率3 GHz、方位角10°附近，由于折射偏轉(zhuǎn)作用，目標底部鏡面反射方向發(fā)生偏轉(zhuǎn)，斜入射下目標RCS增強約10 d B，如圖4（b）與圖4（c）所示。另一方面，隨著頻率的增加，由于等離子體分布噪聲將導(dǎo)致一定的傳播路徑相位起伏，方位角60°附近較大范圍的目標RCS谷值增加，但RCS值仍較小。

由圖4可以看出，在方向角180°的頭部方向，由于鈍錐體結(jié)構(gòu)特點，這個方位本身RCS相對而言較小，加上頭部等離子體較高、吸收作用較大的原因，不同頻段下其RCS減小相對較大。方向角0°和100°方向是兩個明顯的反射方位，其反射較大。方位角100°附近由于入射方位與錐體側(cè)面近似垂直，所以反射較大。當頻率不大于5 GHz的情況，在方位角0°的尾部方向，由于等離子體的折射偏轉(zhuǎn)作用，鈍錐體RCS存在較大減小，甚至完全隱身；隨著頻率增大，等離子體鞘套對目標RCS的影響逐步降低，目標非鏡面反射方向的RCS逐步增加，這是由于等離子體密度的非均勻性導(dǎo)致反射增強。

圖5給出了目標飛行馬赫數(shù)為20時，不同入射電磁波頻率有無等離子體時鈍錐體RCS的對比變化，實線為等離子體覆蓋目標的RCS變化情況，虛線為目標本身的RCS變化情況。等離子體對目標RCS的縮減及增強影響均十分明顯，尤其是入射電磁波頻率為5 GHz時，在較寬方位角范圍內(nèi)等離子體對目標RCS影響明顯。圖5（b）、圖5（c）分別為10 GHz、18 GHz頻率下RCS結(jié)果對比，模擬入射粒子數(shù)增加至60萬。

圖4 飛行馬赫數(shù)為15時不同頻率及方位角的RCS分布

由圖5可以看出，等離子體覆蓋目標后，在方位角20°至80°方向，等離子體鞘套目標RCS較鈍錐體目標本身略大，這主要因為等離子體層局部密度的差異，導(dǎo)致目標光滑表面“粗糙化”，散射相對增強；而在方位角100°至180°方向，由于等離子體吸收作用，等離子體鞘套目標RCS較鈍錐體目標本身RCS縮減。

圖5 飛行馬赫數(shù)為20時鈍錐體不同方位角的RCS分布

4 結(jié)論

針對高超聲速目標等離子體鞘套電磁特性仿真的難點，本文采用場離散的MC方法，仿真得到了等離子體覆蓋目標的電磁特性，并結(jié)合鈍錐體目標及等離子體鞘套目標RCS對比，分析不同頻率下RCS與入射方位角的關(guān)系，分析給出了等離子體對目標RCS的主要縮減及增強機理。MC方法可以對等離子體流場仿真結(jié)果進行自動化處理，便于超高速目標流場-電磁場的一體化仿真建模，可實現(xiàn)等離子體包覆的電大尺寸目標的電磁散射特性的快速計算。