基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法与流程

本发明涉及雷达目标电磁散射特性计算技术领域，具体涉及一种基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法。

背景技术：

常用的等离子体与雷达电磁波作用的仿真方法主要有wkb近似方法、fdtd方法等。wkb方法以几何光学近似为基础，进而求解波动方程，是一种近似方法，只有当等离子体参数随位置的变化足够缓慢时，wkb方法才适用，且一般用于一维非均匀情况。fdtd方法可适用于求解任意分布的飞行器等离子体鞘套，然而当等离子体鞘套分布比较复杂或电磁波频率较高时，fdtd方法需要的计算机内存较大，计算时间较长。而已有蒙特卡罗方法主要用于分析光与等离子体作用的仿真。

代表性参考文献如下：

美国斯坦福大学sriinternational公司的r.j.vidma在1990年利用wkb方法，设计出了一种空间电磁波反射器和吸收器装置(r.j.vidmar,ontheuseofatmospherepressureplasmaaselectro-magneticreflectorandabsorber[j],ieeetransactiononplasmascience,1990,18(4),733-741.)。该研究通过调节等离子体物理参数，实现对不同频率电磁波的吸收和反射，未涉及电磁场粒子化的等离子体电磁特性仿真方法。

美国田纳西州立大学的m.laroussi等人，在1993年发表在ieeetransactiononplasmascience上的文献，给出了基于多层抛物线分布的等离子体密度分布模型(m.laroussi,j.r.roth,numer-icalcalculationofthereflection,absorption,andtransmissionofmicrowavebyanon-uniformplasmaslab[j],ieeetransactiononplasmascience,1993,21(4),366-372.)。该模型采用了wkb方法计算了电磁波在任意角度入射非均匀等离子体情形下的电磁波反射，吸收以及透射情况，未涉及电磁场粒子化的等离子体电磁特性仿真方法。

美国亚利桑那州立大学的j.j.simpson课题组在2015年提出了一种基于随机时域有限差分三维结构电磁波在等离子体中的传播求解(b.t.nguyen,c.furse,j.simpson,a3-dstochasticfdtdmodelofelectromagneticwavepropagationinmagnetizedionos-phereplasma,ieeetransactiononantennapropagation,2015,63(1),304-313.)。该模型由于是采取全波仿真的方法，计算结果虽然准确，但是耗时过长，需要消耗大量的计算机资源，未涉及电磁场粒子化的等离子体电磁特性仿真方法。

美国伊利诺伊州大学香槟分校的k.rajaraman等人在2004年发表在appliedphysics的文章(k.rajaraman,m.kushner,amontecarlosimulationofradiationtrappinginelectrodelessgasdis-chargelamps,1780-1791,vol.37,2004)中采用能量粒子化的蒙特卡罗方法分析了低压等离子体气体灯的发光效率问题，未涉及雷达电磁波与等离子体作用的仿真。

在专利申请方面，西安电子科技大学的石磊等人提出了一种针对动态快速等离子体鞘套电磁波传播计算方法(石磊，姚博，李小平，杨敏，刘彦明，一种快速的动态等离子体鞘套电磁波计算方法，西安电子科技大学，专利申请号：cn201510589009)。其方法中通过建立动态的等离子体鞘套电子密度数学模型，采取了wkb方法获取电磁波的传播特性，避免了直接采取高难度的随机介质建模方法，未涉及电磁场粒子化的等离子体电磁特性仿真方法。

中国人民解放军火箭军工程大学的张金生等人提出了一种磁化等离子体的时域有限差分方法(张金生，李征委，杨东方，一种磁化等离子体的时域有限差分方法，中国人民解放军工程大学，专利号：cn201610157802)。通过已知的等离子体参数模型，针对磁化等离子体，采用把等离子体的物理参数等效为电流源的方法进行电磁波分布的时域有限差分求解，未涉及电磁场粒子化的等离子体电磁特性仿真方法。

中国运载火箭研究院的背景临近空间飞行器系统工程研究所的苏汉生等人提出了一种超高速飞行器等离子体鞘套与电磁波相互作用预测方法(苏汉生，张作一，刘秀祥，李瑾，郑晨，冯数琦，超高速隐身飞行器等离子体鞘套与电磁波相互作用预测方法，北京临近空间飞行器系统工程研究所，专利申请号：cn201611033297)。采用了wkb方法对超高速飞行器再入大气层中的“黑障”问题进行了分析，未涉及电磁场粒子化的等离子体电磁特性仿真方法。

技术实现要素：

本发明的目的在于提供一种基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，解决了三维电大尺寸非均匀等离子体电磁特性难以仿真计算的问题。

一种基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，包含以下步骤：

s1、将照射电磁场粒子化，确定粒子的初始运动趋势；

s2、按某一离散步长，计算该粒子的下一步位置，并从等离子体分布输入模型中提取新位置处等离子体的折射率和吸收系数，其中离散步长要求小于等离子体分布输入模型中最小离散网格尺寸；

s3、根据提取的等离子体等效吸收系数计算该粒子被吸收概率，并用伪随机数判断该粒子是否被吸收，是则跳转至步骤s6，否则执行s4；

s4、根据等离子体折射率的大小及其变化率计算该粒子运动方向的变化；

s5、判断该粒子位置是否到达等离子体模型外廓面，否则跳转至步骤s2开始该粒子的下一步运动循环，是则执行s6；

s6、判断该粒子被吸收或到达等离子体模型外廓面的粒子数是否满足要求，否则跳转至步骤s1开始下一个粒子的仿真循环，是则执行s7；

s7、对达到等离子体外廓面的电磁场粒子进行场积分获取要求的电磁特性。

上述的基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，所述步骤s7中，到达等离子体外廓面的电磁场粒子的散射场公式如下：

其中，λ为电磁波的波长；m为仿真的总粒子数；n为到达离子体外廓面的电磁场粒子数；s表示等离子体输入模型在电磁波入射方向的横截面面积；k为电磁波波数；rtj表示电磁波的发射源点到电磁场粒子起始位置间的距离；rsj表示接收电磁波的观测点到电磁场粒子出射位置间的距离；γj表示电磁场粒子在等离子体中运动的电路径长度，为沿该粒子运动轨迹的等离子体等效折射率的线积分。

上述的基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，所述步骤s6中，要求到达离子体外廓面的电磁场粒子数n与仿真的总粒子数m的比值收敛。

上述的基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，所述步骤s4中，若折射率大于0，电磁场粒子运动方向由下式计算：

折射率不大于0时，电磁场粒子运动方向发生全反射，其运动方向由下式计算：

其中，表示等离子体等效折射率变化的梯度方向；t为在由该粒子上一步运动方向v与决定的平面内垂直于的单位矢量；θ1为矢量v与矢量的夹角。

上述的基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，其中：

δθ≈-tanθ1·δn/n1

其中，n1为该粒子上一步位置处等离子体等效折射率大小；δl为粒子的运动步长。

上述的基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，所述步骤s3中，粒子的吸收概率pa由下式计算：

pa＝1-exp(-a1δl)

其中，a1为该位置处等离子体的等效吸收系数，δl为粒子的运动步长，利用计算机生成0至1间的伪随机数，当所述伪随机数小于pa，则电磁场粒子被吸收。

上述的基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，所述步骤s1中，所述粒子初始运动趋势由下列条件决定：

t1、以输入的等离子体模型外廓面为电磁场粒子运动的起始面；

t2、所述起始面上取随机位置作为粒子起始位置；

t3、所述粒子起始运动方向取电磁波入射方向。

本发明的优点和有益效果是：

(1)直观分析电磁波与非均匀等离子体的作用过程和机理。

(2)提高电磁波与非均匀等离子体的电磁特性的仿真计算效率。

附图说明

图1是本发明中基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法的流程图。

图2是本发明较佳实施例中入射电磁场粒子在等离子体输入模型外廓面分布情况示意图。

图3是本发明较佳实施例中部分电磁场粒子的运动轨迹示意图。

图4是本发明较佳实施例中非均匀等离子体覆盖目标rcs仿真结果。

具体实施方式

以下结合附图，通过详细说明一个较佳的具体实施例，对本发明做进一步阐述。

如图1所示，基于蒙特卡罗方法的非均匀等离子体电磁特性仿真方法，包含以下7个步骤：

s1、将照射电磁场粒子化，确定粒子的初始运动趋势。

粒子初始运动趋势由下列条件决定：

t1、以输入的等离子体模型外廓面为电磁场粒子运动的起始面；

t2、所述起始面上取随机位置作为粒子起始位置；

t3、所述粒子起始运动方向取电磁波入射方向。

如图2所示，本实施例中，等离子体输入模型外廓面尺寸为2.06m×1.63m×0.81m，仿真的总粒子数取m＝50000，粒子起始位置在等离子体输入模型外廓面上随机分布，电磁波取与x轴夹角0°～180°的入射方向。

s2、按某一离散步长，计算该粒子的下一步位置，并从等离子体分布输入模型中提取新位置处等离子体的折射率和吸收系数，其中离散步长要求小于等离子体分布输入模型中最小离散网格尺寸。

取粒子的运动步长δl＝0.01m，计算电磁场粒子位置变化δζ＝vδl，并从等离子体分布参数输入模型中采用线性插值方法提取下一步位置处等离子体的等效折射率和吸收系数。

s3、根据提取的等离子体等效吸收系数计算该粒子被吸收概率，并用伪随机数判断该粒子是否被吸收，是则跳转至步骤s6，否则执行s4。

粒子的吸收概率pa由式(1)计算：

pa＝1-exp(-a1δl)(1)

其中，a1为该位置处等离子体的等效吸收系数，δl为粒子的运动步长，利用计算机生成0至1间的伪随机数，当所述伪随机数小于pa，则电磁场粒子被吸收。

s4、根据等离子体折射率的大小及其变化率计算该粒子运动方向的变化。

若折射率大于0，电磁场粒子运动方向由式(2)计算：

折射率不大于0时，电磁场粒子运动方向发生全反射，其运动方向由式(3)计算：

其中，表示等离子体等效折射率变化的梯度方向；t为在由该粒子上一步运动方向v与决定的平面内垂直于的单位矢量；θ1为矢量v与矢量的夹角。

式(3)中：

δθ≈-tanθ1·δn/n1(4)

其中，n1为该粒子上一步位置处等离子体等效折射率大小；δl为粒子的运动步长。

s5、判断该粒子位置是否到达等离子体模型外廓面，否则跳转至步骤s2开始该粒子的下一步运动循环，是则执行s6。

s6、判断该粒子被吸收或到达等离子体模型外廓面的粒子数是否满足要求，否则跳转至步骤s1开始下一个粒子的仿真循环，是则执行s7。

要求到达离子体外廓面的电磁场粒子数n与仿真的总粒子数m的比值收敛。

仿真的部分粒子运动轨迹如图3所示，其中点划线表示粒子轨迹，离散点表示非均匀等离子体包裹的钝锥体目标外表面。

s7、对达到等离子体外廓面的电磁场粒子进行场积分获取要求的电磁特性。

到达等离子体外廓面的电磁场粒子的散射场公式如下：

其中，λ为电磁波的波长；m为仿真的总粒子数；n为到达离子体外廓面的电磁场粒子数；s表示等离子体输入模型在电磁波入射方向的横截面面积；k为电磁波波数；rtj表示电磁波的发射源点到电磁场粒子起始位置间的距离；rsj表示接收电磁波的观测点到电磁场粒子出射位置间的距离；γj表示电磁场粒子在等离子体中运动的电路径长度，为沿该粒子运动轨迹的等离子体等效折射率的线积分。

按公式(6)，对达到等离子体外廓面的电磁场粒子进行场积分获取钝锥体目标后向rcs随方位角的变化，结果如图4所示，其中方位角0°峰值表示钝锥体底面的镜面反射贡献，方位角100°附近峰值表示钝锥体锥面的镜面反射贡献，方位角180°附近凹谷来源于钝锥体头部较高密度等离子体的吸收衰减。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。