CN104239620B - 一种航天器尾区带电效应仿真方法 - Google Patents

Abstract

本发明提供了一种航天器尾区带电效应仿真方法，构建极区沉降电子能谱的拟合表达式，即极光电子通量Φ随能量E变化的关系，该关系通过功率定律分布Φ_p(E)、麦克斯韦分布Φ_m(E)和高斯分布Φ_G(E)的叠加来表示，并拟合极区沉降电子能谱；建立背景等离子体环境中离子分布特性谱；进行航天器不同表面材料的二次电子发射谱测试；基于上述三类谱采用单元粒子法PIC进行航天器尾区带电效应的过程模拟。具有模拟过程清晰，模拟结果和实际监测带电情况较符合等优点。

Description

一种航天器尾区带电效应仿真方法

技术领域

本发明属于空间计算领域，适用于极轨(Polar Earth Orbit)航天器尾区中表面电介质材料的充电过程模拟，具体涉及一种极区内在背景等离子体和极光电子作用下的航天器尾区带电效应仿真方法。

背景技术

倾角大于或等于55度的低地球轨道(俗称极轨)航天器会频繁穿越极光弧。航天器会遭遇与地球静止轨道相似的高能电子环境。极轨等离子环境主要是增加了由于经过极区从而可能遭遇极光粒子事件，极光粒子事件是指地磁扰动或太阳爆发期间发生的高能带电粒子(电子和质子)沿地磁力线下降到极区引起的极光沉降粒子的增强效应。地球两极处的磁力线由于太阳风的影响，一部分被拉开。这些开力线不再接地磁两极而是有一端通向星际空间，从而在极区形成了一漏斗形状区域，开磁力线延伸到地面极区的部分称之为极盖区。

当卫星运行在低温度、高密度的极区等离子体环境中时，在其尾部形成明显的“航迹”，这是一个不相等的电子和离子耗尽区。由于卫星轨道速率大于离子热速率而小于电子热速率，因此电子可较容易地进入这个区域从而形成一负电位势垒，这就是所谓的“尾迹效应”。它对卫星的明显作用是在尾区介质表面将充电至较高的负电位，此表面电位主要依赖于收集的电子通量和离子通量之比。卫星因尾迹效应而形成的表面不等量带电是影响中低轨道特别是极轨卫星安全运行的重要原因之一。当卫星尾部介质表面带电达到或超过航天器材料击穿阈值后，便会产生静电放电。

目前国际上已有一些在使用的极轨航天器尾区带电效应仿真方法，例如美国NASA的极区大型航天器尾区带电效应仿真方法和ESA的极区航天器尾区等离子体相互作用模拟方法，这些模拟方法原理相似，只是在使用和功能上有些差别，但这些模拟方法的核心代码未公开而且模拟精度无法评估。国内在这方面的研究较少，而且在模拟极区极光电子沉降谱时，未采用较符合实际的极区电子分布谱。

发明内容

有鉴于此，本发明提供了一种航天器尾区带电效应仿真方法，本发明中结合极区沉降电子能谱、离子分布特性谱和材料二次电子发射谱，利用PIC(Particle-In-Cell)方法有效的仿真极区航天器尾区带电效应过程，具有模拟过程清晰，模拟结果和实际监测带电情况较符合等优点。

该航天器尾区带电效应仿真方法，包括如下步骤：

步骤一、构建极区沉降电子能谱的拟合表达式，即极光电子通量Φ随能量E变化的关系，该关系通过功率定律分布Φ_p(E)、麦克斯韦分布Φ_m(E)和高斯分布Φ_G(E)的叠加来表示；利用现有的极区电子谱观测数据获得拟合表达式中的未知参数，得到极区沉降电子能谱表达式；

步骤二、建立背景等离子体环境中离子分布特性谱；

步骤三、进行航天器不同表面材料的二次电子发射谱测试；

步骤四、采用单元粒子法PIC进行航天器尾区带电效应的过程模拟：

利用步骤一中拟合的极区沉降电子能谱，在整个模拟区域所划分的计算网格中分配各能量段的电子个数，然后等效成一定数目的宏粒子；

利用步骤二中建立的离子分布特性谱，在航天器周围空间的计算网格中分配离子密度，然后等效成一定数目的宏粒子；

利用步骤三中获得的每种材料的二次电子发射能谱，分配航天器表面材料附近空间的计算网格中的二次电子密度，然后等效成一定数目的宏粒子；

分配完成后，运行单元粒子法PIC获得模拟结果。

优选地，步骤三对ITO覆膜、Kapton、太阳能电池板玻璃片和碳纤维结构材料进行测试。

优选地，所述步骤二具体为：

步骤1、建立航天器等效模型为表面带负电、半径为a的导体球；a为航天器最大横向半宽；

步骤2、建立离子在势场中轨道运动的能量守恒方程，并对离子质量m和离子电荷q进行无量纲转换，得到无量纲能量守恒方程；

步骤3、建立离子速度空间分布函数f₀()为麦克斯韦速度分布并对速度量进行无量纲转换，则

其中，为单个离子无穷远处的速度矢量，为离子整体运动速度矢量，即航天器速度矢量的反方向；

步骤4、建立离子撞击点处的局部直角坐标系；

局部直角坐标系OXYZ建立在导体球的球面上，原点位于球面的撞击点，Z轴垂直于球面，X在Z平面内且与Z轴垂直，为离子撞击航天器的速度矢量，Y轴满足右手法则；

在局部直角坐标系中定义如下角度：

ψ为航天器速度矢量与在势场中离子撞击航天器的速度矢量之间的夹角，ξ为航天器速度矢量和撞击面法向Z之间的夹角，θ_∞为无穷远处离子的速度与撞击面法向Z之间的夹角，为离子撞击航天器的速度矢量在XY平面上的投影与X轴之间的夹角，即离子运动轨道平面与X轴之间的夹角；

则，上述角度之间的角度关系为：

步骤5、建立球坐标系下的离子充电电流密度表达式：

将局部直角坐标系下的离子充电电流密度表达式转换到球坐标系中，并代入所述无量纲能量守恒方程、所述离子速度空间分布函数和所述角度关系，得到球坐标系下离子充电电流密度表达式：

α和为球坐标系中的极角，对极角积分得到：

其中，

x＝-v₀V₀sinξsinθ_∞ (10)

其中，v₀为的大小，V₀为的大小；

步骤6、求解贝塞尔函数从而消去充电电流密度表达式(6)中的未知量；为

步骤7、根据航天器等效模型，建立势场中离子的运动轨道方程，并利用边界条件获得角度θ₀、θ_∞和α的关系为：

将(12)式代入(13)式得到θ_∞与α的关系式I；

步骤8、将θ_∞与α的关系式I代入式(7)，从而消去θ_∞，得到j与ξ、α的关系式II；接着利用关系式II对α进行积分，得到j与ξ的关系式III；

步骤9、对j与ξ的关系式III中的无量纲化速度进行有量纲化转换，得到即航天器离子分布特性；

则步骤四中，所述利用步骤二中建立的离子分布特性谱，在航天器周围空间的计算网格中分配离子密度具体为：根据航天器离子分布特性，分配航天器周围空间网格各个角度上的离子充电电流密度。

有益效果：

(1)本发明中极区沉降电子谱可通过三种分布的叠加来表示，实现极光电子谱的拟合，拟合结果符合美国DMSP卫星实测谱。

(2)本发明实施例背景等离子体环境中离子分布特性谱计算中，把航天器表面材料等效成同种材料的并处于一定负电位的球面，利用轨道限制理论计算出的离子电流体现出随撞击点角度、离子整体速度和球表面电势的变化，结果符合实际充电过程。

附图说明

图1为极光电子沉降谱拟合结果。

其中，点线-功率定律分布谱，点虚线-麦克斯韦分布谱，虚线-高斯分布谱，实线-极光电子拟合谱。

图2为航天器尾区带电效应仿真结果。

图3(a)和图3(b)为局部坐标系示意图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提供了一种尾区带电效应仿真分析方法，结合极区沉降电子能谱、离子分布特性谱和材料二次电子发射谱，利用PIC方法能够有效的仿真分析极区航天器尾区带电效应过程。该方案具体包括下列步骤：

步骤一、通过三种分布的叠加来表示极区沉降电子能谱。

构建极区沉降电子能谱的拟合表达式，即极光电子通量Φ随能量E变化的关系，该关系通过功率定律分布Φ_p(E)、麦克斯韦分布Φ_m(E)和高斯分布Φ_G(E)的叠加来表示。其中，功率定律分布函数为Φ_p＝A_pE^-α1，式中A_p为系数，α₁为指数，E为能量。麦克斯韦分布函数为式中e为单位电子电荷，m_e为电子的质量，n为电子密度，θ₁为能量参数，E为能量。高斯分布为式中A_G为系数，E₀和Δ为能量参数，E为能量。

利用现有的极区电子谱观测数据获得拟合表达式中的未知参数，得到极区沉降电子能谱表达式。

本发明实施例中功率定律分布与二次电子和背散射电子的原始沉降电子束能量相关，麦克斯韦分布能谱代表具有中、高能(用温度表示)的环境束流，高斯分布代表极光增强的电子流。对于大多数环境下，上述分布谱中的参数设置如下：A_p＝3×10¹¹ _m ^-3，α₁＝1.1，E_PL＝50eV，E_PH＝1.6×10⁶eV，n＝6×10⁵m^-3，θ₁＝8keV，A_G＝4×10₄m^-3，E₀＝24keV，Δ＝16keV。

步骤二、建立背景等离子体环境中离子分布特性谱。

可以按照现有技术，利用麦克斯韦-玻尔兹曼分布来表达背景等离子体环境中离子分布特性谱。

步骤三、进行航天器不同表面材料的二次电子发射谱测试。

本实施例中分别选用航天器表面材料Kapton、ITO、太阳能电池板玻璃片和材料在二次电子发射系数测试设备中测试，分别获得其二次电子发射能谱，为PIC数值模拟建立基础。

步骤四、PIC(单元粒子法)方法过程模拟；

利用步骤一中拟合的极区沉降电子能谱，在图2中整个模拟区域所划分的计算网格中分配各能量段的电子个数，然后等效成一定数目的宏粒子；

利用步骤二中建立的离子分布特性谱，在航天器周围空间的计算网格中分配离子密度，然后等效成一定数目的宏粒子；

利用步骤三中获得的每种材料的二次电子发射能谱，分配航天器表面材料附近空间的计算网格中的二次电子密度，然后等效成一定数目的宏粒子；

航天器模型采用在轨的极轨对地观测卫星缩比模型，其显著特点是拥有较大尺寸的抛物面天线。航天器表面材料设置采用数据库方式直接调用。航天器周围空间模拟区域网格划分采用均匀网格划分可以简化程序算法；

基于上述设置，运行单元粒子法PIC获得模拟结果。PIC的具体步骤：获取离子和电子的位置和速度后，将它们所带的电荷和电流分配到网格节点上，求得网格节点处的电荷密度和电流密度，然后求解Maxwell方程即可得到网格点上电磁场值，再通过插值得到离子和电子处的电磁场，粒子在Lorentz力的作用下运动。通过如此的循环迭代过程自洽求解等离子体物理问题。

上述实施例中的步骤二中，采用常规的麦克斯韦-玻尔兹曼分布来表达背景等离子体环境中离子分布特性谱。但实际上，麦克斯韦-玻尔兹曼分布对极区内背景离子分布特性描述中，一般认为航天器周围离子分布都是均匀的，不考虑离子充电电流在航天器撞击面和尾区面的变化。然而实际上，航天器尾区不同位置处或者说不同角度上，充电电流是不同的，因此有必要建立中低轨航天器表面充电电流密度与某个能够表达离子入射位置的参量的计算模型，利用该模型来模拟航天器尾区充电过程，从而提高模拟和仿真精度。

为此，本发明进一步提供了一种航天器离子分布特性的计算方法，该方法利用轨道限制理论近似计算离子的充电电流，即在等效计算球静止坐标系中半径为a、表面处于负电位的导体球表面离子的充电电流，计算出的离子电流体现出随撞击点角度、离子整体速度和球表面电势的变化，结果符合实际充电过程。

该方案具体包括下列步骤：

步骤1、建立航天器等效模型：

设航天器最大横向半宽为a，建立航天器等效模型为表面带负电、半径为a的导体球。

航天器等效模型的构建原理为：

①轨道环境中离子撞击航天器表面的过程实际上是带电粒子在势场中的运动过程，因此可以将航天器认为是一个带负电的物体；

②离子撞击航天器表面后被航天器吸收，形成充电电流，而不是反弹出去，所以将航天器认为是一个导体；

③本发明希望得到离子撞击航天器不同位置时产生的充电电流密度j，而采用球形模拟航天器，可以用撞击点法向与航天器运动方向的夹角来表示不同的撞击点位置，这样容易计算。

基于上述分析，本实施例将航天器简化为带负电的导体球。

步骤2、建立离子在势场中轨道运动的能量守恒方程。

能量守恒方程为：

其中，m为离子质量，v为单个离子速度，q为离子电荷，Φ为导体球表面电势，E为离子总能量，v₀为单个离子无穷远处的速度。

本实施例为了简化计算过程，将能量守恒方程进行无量纲转换，即令q＝m＝1，从而得到无量纲能量守恒方程为：

后面步骤得到的j与ξ的关系式III中没有出现m和q，因此不需要进行后续的有量纲化处理。

步骤3、建立不考虑势场情况下的离子速度空间分布函数f₀()。

如果实际环境中没有外界干扰，即不考虑势场情况下物体的相互势场影响，则离子处于热平衡状态，麦克斯韦分布描述了这种状态，因此离子速度空间分布符合麦克斯韦分布，所建立的离子速度空间分布函数f₀()为：

其中

为了简化计算，对上式中的速度量进行无量纲化处理，即去除v_T，从而上式转换为：

其中，为单个离子无穷远处的速度矢量，为离子整体运动速度矢量(即航天器运行速度矢量的反方向)，v_T为单个离子的热运动速度，k为玻尔兹曼常量，T为离子温度。

步骤4、建立撞击点局部直角坐标系；

局部直角坐标系OXYZ建立在导体球的球面上，原点位于球面的撞击点，Z轴垂直于球面，X在平面内且与Z轴垂直，为势场中单个离子的速度矢量，在该局部直角坐标系中就是离子撞击航天器的速度矢量，Y轴满足右手法则。

在局部直角坐标系中定义如下参数：

ψ为航天器速度()与在势场中离子撞击航天器的速度矢量之间的夹角，ξ为航天器速度()和撞击面法向Z之间的夹角，θ_∞为无穷远处离子的速度与撞击面法向Z之间的夹角，为离子撞击航天器的速度矢量在XY平面上的投影与X轴之间的夹角，即离子运动轨道平面与X轴之间的夹角。

基于上述角度定义，如图1，角度关系满足如下等式：

步骤5、在局部直角坐标系下，导体球表面撞击点处法向方向的离子充电电流密度为：

其中，为势场中单个离子的速度矢量，为撞击点外法向方向的单位矢量。

将离子充电电流密度表达式(5)从局部直角坐标系变换到球坐标系，球坐标系中的参量为极角α和，α∈[0,90°]，∈[0,360°]，并且将式(2)(3)(4)代入(5)，则球坐标系下离子充电电流密度为：

对极角积分得到：

其中，

x＝-v₀V₀sinξsinθ_∞ (10)

其中，v₀为的大小，V₀为的大小；

式(7)中包括4个未知量：α、、ξ和θ_∞，下面通过数值求解的方式消去α、和θ_∞，最终得到j与ξ的关系式，即离子分布特性。但是由于j的表达式过于复杂，因此本发明涉及了如下步骤六～八来简化求解过程，从而降低求解难度。

步骤6、恰好为零级贝塞尔函数的表达形式，可以利用幂级数展开实现数值求解，因此本步骤定义中的xcos为x'，对x'进行幂级数展开，取前7项分别进行的积分，从而消去了。

中有指数积分部分，其计算过程太过复杂，因此本发明将其转换为求和积分，从而简化积分计算。经本步骤的处理，四个未知量仅剩下三个。

步骤7、根据航天器等效模型，建立势场中离子的运动轨道方程为：

其中，r为离子轨道半径，θ₀为轨道参数。

利用边界条件可获得角度θ₀、θ_∞和α的关系：

1)当r＝a时，θ＝0，运动轨道方程变为：

该式表达出θ₀与α的关系；

2)当r→∞时，θ＝θ_∞，运动轨道方程变为：

将(12)式代入(13)得到θ_∞与α的关系式I。

步骤8、将θ_∞与α的关系式I代入式(7)，从而消去θ_∞，得到j与ξ、α的关系式II；接着利用关系式II对α进行0～π/2的积分，得到j与ξ的关系式III。图3绘出了关系式III的曲线，横轴为撞击点角度ξ，纵轴为j。

步骤9、对j与ξ的关系式III进行有量纲化转换，得到即航天器离子分布特性。航天器离子分布特性表达了有量纲的j与角度ξ的关系。

由于步骤3中对速度进行了无量纲化处理，且最终得到关系式III中存在无量纲化的速度，因此需要将其进行有量纲化处理，即将j乘以v_T，即为离子通量，最终得到航天器离子分布特性。

至此，本流程结束。

采用上述步骤9获得的航天器离子分布特性，在航天器周围空间的计算网格中分配离子密度具体为：根据航天器离子分布特性，分配航天器周围空间网格各个角度上的离子充电电流密度。这样一来，仿真结果更为真实和准确。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种航天器尾区带电效应仿真方法，其特征在于，包括如下步骤：

步骤一、构建极区沉降电子能谱的拟合表达式，即极光电子通量Φ随能量E变化的关系，该关系通过功率定律分布Φ_p(E)、麦克斯韦分布Φ_m(E)和高斯分布Φ_G(E)的叠加来表示；利用现有的极区电子谱观测数据获得拟合表达式中的未知参数，得到极区沉降电子能谱表达式；

步骤二、建立背景等离子体环境中离子分布特性谱；

步骤三、进行航天器不同表面材料的二次电子发射谱测试；

步骤四、采用单元粒子法PIC进行航天器尾区带电效应的过程模拟：

利用步骤一中拟合的极区沉降电子能谱，在整个模拟区域所划分的计算网格中分配各能量段的电子个数，然后等效成一定数目的宏粒子；

利用步骤二中建立的离子分布特性谱，在航天器周围空间的计算网格中分配离子密度，然后等效成一定数目的宏粒子；

利用步骤三中获得的每种材料的二次电子发射能谱，分配航天器表面材料附近空间的计算网格中的二次电子密度，然后等效成一定数目的宏粒子；

分配完成后，运行单元粒子法PIC获得模拟结果；

所述步骤二具体为：

步骤1、建立航天器等效模型为表面带负电、半径为a的导体球；a为航天器最大横向半宽；

步骤2、建立离子在势场中轨道运动的能量守恒方程，并对离子质量m和离子电荷q进行无量纲转换，得到无量纲能量守恒方程；

步骤3、建立离子速度空间分布函数为麦克斯韦速度分布并对速度量进行无量纲转换，则

其中，为单个离子无穷远处的速度矢量，为离子整体运动速度矢量，即航天器速度矢量的反方向；

步骤4、建立离子撞击点处的局部直角坐标系；

局部直角坐标系OXYZ建立在导体球的球面上，原点位于球面的撞击点，Z轴垂直于球面，X在Z平面内且与Z轴垂直，为离子撞击航天器的速度矢量，Y轴满足右手法则；

在局部直角坐标系中定义如下角度：

ψ为航天器速度矢量与在势场中离子撞击航天器的速度矢量之间的夹角，ξ为航天器速度矢量和撞击面法向Z之间的夹角，θ_∞为无穷远处离子的速度与撞击面法向Z之间的夹角，为离子撞击航天器的速度矢量在XY平面上的投影与X轴之间的夹角，即离子运动轨道平面与X轴之间的夹角；

则，上述角度之间的角度关系为：

步骤5、建立球坐标系下的离子充电电流密度表达式：

将局部直角坐标系下的离子充电电流密度表达式转换到球坐标系中，并代入所述无量纲能量守恒方程、所述离子速度空间分布函数和所述角度关系，得到球坐标系下离子充电电流密度表达式：

α和为球坐标系中的极角，对极角积分得到：

$j={\left(2\mathrm{\π}\right)}^{-3/2}{\∫}_0^{\mathrm{\∞}}{v}_0(v_0^2+2\mathrm{\Φ})\exp \[-v_0^2/2\]F\left(v_0\right){\mathrm{dv}}_0---\left(7\right)$

其中，

$F\left(v_0\right)={\left(2\mathrm{\π}\right)}^{-3/2}{\∫}_0^{\mathrm{\π}/2}\exp \[-V_0^2/2\]\exp \[\left(v_0{V}_0{\mathrm{cos\ξcos\θ}}_{\mathrm{\∞}}\right)/2\]I_0\left(x\right)cos\mathrm{\α}\sin \mathrm{\α}d\mathrm{\α}---\left(8\right)$

x＝-v₀V₀sinξsinθ_∞ (10)

其中，v₀为的大小，V₀为的大小；

步骤6、求解贝塞尔函数从而消去充电电流密度表达式(6)中的未知量

步骤7、根据航天器等效模型，建立势场中离子的运动轨道方程，并利用边界条件获得角度θ₀、θ_∞和α的关系为：

${\mathrm{cos\θ}}_0=\[2(E/\mathrm{\Φ}+1){\sin }^2\mathrm{\α}-1\]{\sqrt{1+4E(E+\mathrm{\Φ})/{\mathrm{\Φ}}^2\·{\sin }^2\mathrm{\α}}}^{-1}---\left(12\right)$

$cos({\mathrm{\θ}}_{\mathrm{\∞}}-{\mathrm{\θ}}_0)=-{\sqrt{1+4E(E+\mathrm{\Φ})/{\mathrm{\Φ}}^2\·{\sin }^2\mathrm{\α}}}^{-1}---\left(13\right)$

将(12)式代入(13)式得到θ_∞与α的关系式I；

步骤8、将θ_∞与α的关系式I代入式(7)，从而消去θ_∞，得到j与ξ、α的关系式II；接着利用关系式II对α进行积分，得到j与ξ的关系式III；

步骤9、对j与ξ的关系式III中的无量纲化速度进行有量纲化转换，得到即航天器离子分布特性；

则步骤四中，所述利用步骤二中建立的离子分布特性谱，在航天器周围空间的计算网格中分配离子密度具体为：根据航天器离子分布特性，分配航天器周围空间网格各个角度上的离子充电电流密度。

2.如权利要求1所述的方法，其特征在于，步骤三对ITO覆膜、Kapton、太阳能电池板玻璃片和碳纤维结构材料进行测试。