CN109979543A - 用于高密度、大尺寸等离子体的粒子模拟方法 - Google Patents

Abstract

一种用于高密度、大尺寸等离子体的粒子模拟方法，首先建立电推进推力器等离子体的粒子模拟初始模型，选取适当的缩比因子。分别针对无碰撞粒子运动、等离子体粒子碰撞过程，推导缩比前后模型中涉及到的物理量变化关系，得到缩比模型中的物理量。执行粒子模拟方法，计算得到缩比模型中的等离子体流动分布以及相关结果。最后，通过本发明中的方法逆推得到原始状态等离子体流动特性相关结果。本发明能够实现高密度、大尺寸下的等离子体粒子模拟，且能够有效减小总计算量、减小模拟中达到平衡时间，提高等离子体粒子模拟方法的运算效率。

Description

用于高密度、大尺寸等离子体的粒子模拟方法

技术领域

本发明属于等离子体的数值模拟技术领域，涉及等离子体的粒子模拟方法，具体来说，是一种用于高密度、大尺寸等离子体的粒子模拟方法。

背景技术

粒子模拟方法利用计算机跟踪大量单个粒子的运动，是模拟等离子体运动的重要方法之一。由于这种模拟方法是基于等离子体中大量真实粒子的运动而非将等离子体视为连续的流体介质，因此，其准确性远高于其它方法，并且能够反映等离子体的微观特性，可以用粒子动力学方式精确描述等离子体。粒子模拟方法依靠其在等离子体模拟中的优势，广泛应用于天体射流、磁层物理、热核聚变以及电推进等研究领域。但是，由于粒子模拟方法以单个粒子为模拟单元的特点，模拟过程中对时间步长Δt、网格长度Δx等存在较为苛刻的限制：

式中，c代表光速，ω_pe代表等离子体电子频率，λ_De代表等离子体德拜长度。以电推进磁等离子体推力器为例，根据推力器等离子体的粒子数密度(高达10²¹/m³)与电子温度等实验测量数据估计，等离子体电子振荡频率和德拜长度的数量级分别为10¹²Hz和10^-6m。一般情况下，推力器的放电通道长度量级为10^-2m，显然每个方向至少需要大约10⁴个网格才能满足上述限制条件，如果为二维粒子模拟模型，则超过10⁸个计算网格。另一方面，模拟时间取决于中性粒子(速度量级10²m/s)通过放电通道的运动时间，时间量级为10^-4s。因此，在时间步长量级为10^-12s的情况下，需要超过10⁸个运算时间步才能达到平衡状态。考虑到上述高密度、大尺寸下的等离子体粒子模拟计算量如此之大，除非采取合理的计算加速方法，否则采用粒子模拟方法是不切实际的。

为了实现高密度、大尺寸下的等离子体粒子模拟，学者们提出了许多不同的方法。比如，2001年Szabo提出两种途径，通过减小重粒子的质量(文献Fully kinetic numericalmodeling of a plasma thruster.Szabo,PhD.thesis,Massachusetts Institute ofTechnology,Boston,2001)，相应地增加速度，进而可以减少重粒子在放电通道中的运动时间，使达到平衡的计算时间缩短；通过提出人为介电常数代替计算中的真空介电常数，增加德拜长度和时间步长，进而放宽对时间步长与网格长度的限制。该加速方法能够从一定程度上提高计算效率，但是在实际应用中存在以下缺点：

1、减小重粒子质量，会改变原来的离子电流与电子电流比例，粒子模拟结果不能反映出等离子体的原始状态；

2、改变真空介电常数并不能改变磁场约束下电子的回旋频率，时间步长仍然足够小，导致实际的加速效果并不明显。

针对上述方法的不足，在2005年，Taccogna提出在保证推力或比冲等重要参数不变的情况下通过等比例缩小计算模型的尺寸来减小计算量(文献Self-similarity inHall plasma discharges:Applications to particle models.Taccogna,physics ofplasmas,2005)。然而，该方法同样带来一些难以解决的问题，首先，该方法不能从粒子轨迹以及碰撞效果方面还原原始放电状态，同时，由于没有适当修正缩比模型中粒子能量与碰撞截面的关系，导致该方法对等离子体放电参数有一定的范围要求，比如，电压太小的情况下缩比模型中电子无法正常电离气体导致不能模拟放电过程，限制了这种方法的应用。

发明内容

针对现有技术无法高效、精确地对电推进推力器中密度较高、尺度较大的放电等离子体进行粒子模拟的问题，本发明提供一种用于高密度、大尺寸等离子体的粒子模拟方法。该方法针对高密度、大尺寸等离子体粒子模拟中存在巨大计算量的难题，提出一种缩比模型的思想，将推力器几何尺寸等比例缩小，在此过程中保证单粒子的运动轨迹在模型缩小前后相对不变，保证粒子碰撞过程在模型缩小前后保持不变。本发明能够实现高密度、大尺寸下等离子体的粒子模拟，并且能够有效减小总计算量、减小模拟中达到平衡时间，提高等离子体粒子模拟方法的运算效率。

为实现上述技术目的，本发明采用的技术方案如下：

用于高密度、大尺寸等离子体的粒子模拟方法，包括以下步骤：

第一步，建立推力器等离子体的粒子模拟初始模型；

第二步，选取缩比因子，将推力器等离子体的粒子模拟初始模型中的推力器通道尺寸按照缩比因子等比例缩小，缩小后对应的模拟模型即推力器等离子体的粒子模拟的缩比模型；

第三步，根据带电粒子在电磁场中的运动轨迹，确定缩比前后模型中无碰撞粒子运动状态物理量的缩比关系；

第四步，保持缩比前后模型中粒子碰撞过程不变，确定缩比前后模型中等离子体粒子碰撞过程物理量的缩比关系；

第五步，通过第三步和第四步中得到的缩比关系，确定缩比前后模型中其它物理量的缩比关系；

第六步，根据第三步、第四步和第五步得到的各物理量的缩比关系，得到缩比模型中的各物理量；执行粒子模拟中电磁场求解以及粒子运动过程求解，计算得到缩比模型中等离子体的运动状态，得到等离子体的流动特性。

第七步：根据第六步中得到的计算结果以及第三步、第四步和第五步得到的各物理量的缩比关系，逆推得到推力器初始模型中等离子体的运动状态以及流动特性。

本发明第一步中，推力器等离子体的粒子模拟初始模型如下：

式中，E、B分别为电场与磁场，J为电流密度，ρ_q为电荷密度，ε₀为真空介电常数，μ₀为真空磁导率，m、V、q分别为粒子质量、速度与电荷量。

本发明第二步中，缩比因子ζ∈(0，1)。

本发明第三步中，包括确定缩比前后模型中的带电粒子在磁场中的粒子回旋半径的缩比关系、缩比前后模型中电场、磁场和粒子速度的缩比关系等，方法如下：

缩比因子为ζ，则有：

L_scal＝ζL (2)

式中，L为推力器通道尺寸，L_scal为缩比模型中推力器通道尺寸；

为了能够保证粒子相对运动轨迹在缩比前后模型中的不变性，做出以下规定：

缩比前后磁场对带电粒子的约束效果不变：

式中，r为粒子回旋半径；

为保证粒子轨迹相对位置的固定，缩比前后模型中的带电粒子在磁场中的粒子回旋半径的缩比关系与缩比前后模型对应的推力器通道尺寸的缩比关系相同，即：

r～L (4)

因此，得到：

缩比前后模型中带电粒子在电磁场中的漂移轨迹相对位置保持不变，即认为漂移运动受到一个向心力作用并保持不变：

其中，F_cyl为向心力，v_DE为粒子漂移速度，R为圆周运动半径；

得到：

粒子的能量完全来自电场做功：

1/2mv²＝qE·R (8)

得到：

υ～(EL)^1/2 (9)

联立式(5)(7)(9)得到缩比前后模型中电场、磁场和粒子速度的缩比关系为：

本发明第四步中，包括确定缩比前后模型中碰撞截面σ、粒子动能E_th以及粒子数密度n的缩比关系，方法如下：

维持缩比前后模型中放电特征不改变，则必须要保证克努森数与霍尔参数在缩比前后模型中保持一致即均为一个常数：

式中，const表示常数，λ_ij为粒子碰撞自由程，ω_c为电子回旋频率，υ_ij为等离子体中i与j粒子之间碰撞频率；

进而得到：

式中，n为粒子数密度，σ为等离子体中i与j粒子之间的碰撞截面；

缩比前后模型中碰撞截面σ保持不变，缩比前后模型中粒子动能E_th以及粒子数密度n做相应变化，缩比前后模型中碰撞截面σ、粒子动能E_th以及粒子数密度n的缩比关系为：

式中，E_th为缩比模型中碰撞截面对应的粒子动能。

本发明第五步中，包括确定缩比前后模型中推力器放电电流的缩比关系以及缩比前后模型中推力器推力和推进效率的缩比关系，具体如下：

缩比前后模型中推力器放电电流的变化关系为：

I_i，e～L^3/2

缩比前后模型中推力器推力和推进效率的变化关系为：

F～L²，η～1。

利用本发明可以在保证粒子运动状态、碰撞效果完全不变的情况下，得到本发明中缩比模型的物理量变化关系。相对于原始粒子模拟过程，采用本发明提出的缩比模型模拟的有益效果主要体现在以下几方面：

(1)减少模拟的总粒子数。假设二维模拟中某种放电状态下达到平衡时粒子的平均密度为N_p，对应的模型尺寸为L，则模拟的粒子总数量为N＝N_pL²。缩比模型中模拟的粒子总数为N_ζ＝ζN_pL²，即N_ζ＝ζN。类似地，如果三维模型模拟，则模拟的粒子总数量与原始粒子数关系为N_ζ＝ζ²N。

(2)减少达到平衡的时间。假设某电子从阴极发处经历时间t后跑出模拟区域(被模拟边界吸收)，由于粒子运行时间变化规律为则在缩比模型中具有相同路径的电子经历的模拟时间为t_ζ＝ζ^1/2t，从而使得模拟可以更快的达到平衡状态。

(3)能够还原放电状态。由于本发明方法的物理量变化关系推导过程完全保证了粒子运动轨迹以及粒子碰撞过程在模型缩比前后保持不变，因此，本发明方法是自洽的，模拟结果可以反应出等离子体分布、粒子组成以及碰撞过程相关的全部原始放电状态。

(4)本发明方法的适用范围广泛。适用于不同情况外加电磁场中的等离子体粒子模拟过程，包括复杂电磁场作用下，也可以用于单一电场或磁场作用下的等离子体模拟。

附图说明

图1为本发明的流程图。

图2展示的是一实施例中采用本发明方法计算得到的电子在推力器通道中的分布情况图。

具体实施方式

下面结合附图对本发明的实施方式作详细说明。

参照图1，为本发明一实施例的流程图。一种用于高密度、大尺寸等离子体的粒子模拟方法，首先建立推力器等离子体的粒子模拟初始模型，选取适当的缩比因子。分别针对无碰撞粒子运动、等离子体粒子碰撞过程，推导缩比前后模型中涉及到的物理量变化关系，得到缩比模型中的物理量。执行粒子模拟方法，计算得到缩比模型中的等离子体流动分布以及相关结果。最后，通过本发明中的方法逆推得到原始状态等离子体流动特性相关结果。具体步骤如下：

第一步，建立推力器等离子体的粒子模拟初始模型。

为了简化计算，模型中忽略麦克斯韦方程中的位移电流，电磁场以及粒子运动模型为：

式中，E、B分别为电场与磁场，J为电流密度，ρ_q为电荷密度，ε₀为真空介电常数，μ₀为真空磁导率，m、V、q分别为粒子质量、速度与电荷量。

第二步，确定缩比因子ζ∈(0，1)。

根据经验结果或已有实验结果，计算等离子体频率以及德拜长度，可由下式得出：

式中，e、m_e、n_e、T_e分别为电子电荷量、质量、数密度与温度。

根据上述结果，计算时间步长Δt以及空间网格长度Δx，其中必须满足粒子模拟中关于时间空间尺度的要求，如式(3)：

c代表光速，ω_pe代表等离子体电子频率，λ_De代表等离子体德拜长度。

推力器等离子体的粒子模拟初始模型缩比前对应的推力器通道尺寸L为已知条件，根据推力器通道尺寸L与计算得到的空间网格长度Δx，确定缩比因子ζ∈(0，1)。

第三步，根据带电粒子在电磁场中的运动轨迹，确定缩比前后模型中无碰撞粒子运动状态物理量的缩比关系。

缩比因子为ζ∈|0，1]，则有：

L_scal＝ζL (2)

式中，L为推力器通道尺寸，L_scal为缩比模型中推力器通道尺寸。

为了能够保证粒子相对运动轨迹在缩比前后模型中的不变性，做出以下规定：

缩比前后磁场对带电粒子的约束效果不变：

式中，r为粒子回旋半径。

为保证粒子轨迹相对位置的固定，缩比前后模型中的带电粒子在磁场中的粒子回旋半径的缩比关系与缩比前后模型对应的推力器通道尺寸的缩比关系相同，即：

r～L (4)

因此，可以得到：

缩比前后模型中带电粒子在电磁场中的漂移轨迹相对位置保持不变，即可以认为漂移运动受到一个向心力作用并保持不变：

其中，F_cyl为向心力，v_DE为粒子漂移速度，R为圆周运动半径。

得到：

粒子的能量完全来自电场做功：

1/2mv²＝qE·R (8)

得到：

υ～(EL)^1/2 (9)

联立式(5)(7)(9)得到缩比后的模型中电场、磁场和粒子速度的缩比关系：

第四步，保持缩比前后模型中粒子碰撞过程不变，确定缩比前后模型中等离子体粒子碰撞过程物理量的缩比关系。

克努森数表征等离子体的稀薄程度，其值为粒子平均自由程与放电腔特征长度的比值。霍尔参数指的是在磁场作用下，带电粒子回旋频率与碰撞频率的比值。如果要维持缩比前后模型的放电特征不改变，则必须要保证克努森数与霍尔参数在缩比前后模型中保持一致：

式中，λ_ij为粒子碰撞自由程，ω_c为电子回旋频率，υ_ij为等离子体中i与j粒子之间碰撞频率。

进而得到：

式中，n为粒子数密度，σ为等离子体中i与j粒子之间的碰撞截面。

由于碰撞截面σ与粒子的动能相关，在第三步中已知推力器缩比后模型中粒子速度的变化，因此可以认为模型缩比前后碰撞截面保持不变，粒子速度(动能)在缩比后模型中做出相应变化，结合式(10)可得到，缩比后的模型中，碰撞截面σ、粒子动能E_th以及粒子数密度n与推力器通道尺寸变化关系为：

式中，E_th为缩比模型中碰撞截面对应的粒子动能。

第五步，通过第三步和第四步中得到的缩比关系，确定缩比前后模型中其它物理量的缩比关系；

推力器放电电流：

I_i，e＝en_i，eυ_i，eS (14)

式中S为放电通道横截面积，v_i,e为带电粒子速度，n_i,e为离子与电子数密度；

显然，缩比后模型中推力器放电电流的变化关系为：

I_i，e～L^3/2 (15)

推力器推力、推力器的推进效率：

对稳态工作的推力器，由动量定理可得：

式中F_h为推力器合外力，m为推力系统总质量，假设真空中不受合外力作用，即F_h＝0，为推力器的推力，则有：

式中，η为推力器的推进效率，F为推力器推力，v代表离子喷出推力器放电腔的速度，I和V分别为推力器的工作电流和电压。

由式(17)可知，缩比后模型中推力器的推力和推进效率的变化关系为：

F～L²，η～1 (18)

第六步，根据第三步、第四步和第五步得到的各物理量的缩比关系，得到缩比模型中的各物理量；执行粒子模拟方法中电磁场求解以及粒子运动过程求解，计算得到缩比模型中等离子体的运动状态，得到等离子体的流动特性。可参见：等离子体粒子模拟，邵福球，科学出版社，2002。

第七步：根据第六步中得到的计算结果以及第三步、第四步和第五步得到的各物理量的缩比关系，逆推得到原始推力器等离子体的粒子模拟仿真模型中等离子体的运动状态以及流动特性。图2展示的是采用本发明方法计算得到的电子在推力器通道中的分布。

本发明中的方法通过采用模型等比例缩小的方式，解决了密度较高、尺度较大的放电等离子体粒子模拟中存在的计算量大、计算耗时长等难题。由于本发明方法是在基于单粒子的运动状态以及粒子碰撞过程在模型缩小前后相对不变的情况下得到的，因此本发明方法是自洽的，模拟结果可以完全反应出原始等离子体的状态。本发明方法能够有效减少总模拟粒子数、减小单粒子在模型中的运动时间，从而在减少总计算量以及模拟达到平衡时间两个方面起到加速计算的效果，提高等离子体粒子模拟方法的运算效率。

以上所述仅为本发明的优选的实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.用于高密度、大尺寸等离子体的粒子模拟方法，其特征在于：包括以下步骤：

第一步，建立推力器等离子体的粒子模拟初始模型；

第二步，选取缩比因子，将推力器等离子体的粒子模拟初始模型中的推力器通道尺寸按照缩比因子等比例缩小，缩小后对应的模拟模型即推力器等离子体的粒子模拟的缩比模型；

第三步，根据带电粒子在电磁场中的运动轨迹，确定缩比前后模型中无碰撞粒子运动状态物理量的缩比关系；

第四步，保持缩比前后模型中粒子碰撞过程不变，确定缩比前后模型中等离子体粒子碰撞过程物理量的缩比关系；

第五步，通过第三步和第四步中得到的缩比关系，确定缩比前后模型中其它物理量的缩比关系；

第六步，根据第三步、第四步和第五步得到的各物理量的缩比关系，得到缩比模型中的各物理量；执行粒子模拟中电磁场求解以及粒子运动过程求解，计算得到缩比模型中等离子体的运动状态，得到等离子体的流场特性。

第七步：根据第六步中得到的计算结果以及第三步、第四步和第五步得到的各物理量的缩比关系，逆推得到原始推力器等离子体的粒子模拟仿真模型中等离子体的运动状态以及流动特性。

2.根据权利要求1所述的用于高密度、大尺寸等离子体的粒子模拟方法，其特征在于：第一步中，推力器等离子体的粒子模拟初始模型如下：

式中，E、B分别为电场与磁场，J为电流密度，ρ_q为电荷密度，ε₀为真空介电常数，μ₀为真空磁导率，m、V、q分别为粒子质量、速度与电荷量。

3.根据权利要求1或2所述的用于高密度、大尺寸等离子体的粒子模拟方法，其特征在于：第二步中，缩比因子ζ∈(0，1)。

4.根据权利要求3所述的用于高密度、大尺寸等离子体的粒子模拟方法，其特征在于，第三步中：

缩比因子为ζ，则有：

L_scal＝ζL(2)

式中，L为推力器通道尺寸，L_scal为缩比模型中推力器通道尺寸；

为了能够保证粒子相对运动轨迹在缩比前后模型中的不变性，做出以下规定：

缩比前后磁场对带电粒子的约束效果不变：

式中，r为粒子回旋半径；

为保证粒子轨迹相对位置的固定，缩比前后模型中的带电粒子在磁场中的粒子回旋半径的缩比关系与缩比前后模型对应的推力器通道尺寸的缩比关系相同，即：

r～L(4)

因此，得到：

缩比前后模型中带电粒子在电磁场中的漂移轨迹相对位置保持不变，即认为漂移运动受到一个向心力作用并保持不变：

其中，F_cyl为向心力，v_DE为粒子漂移速度，R为圆周运动半径；

得到：

粒子的能量完全来自电场做功：

1/2mv²＝qE·R(8)

得到：

v～(EL)^1/2(9)

联立式(5)(7)(9)得到缩比前后模型中电场、磁场和粒子速度的缩比关系为：

5.根据权利要求4所述的用于高密度、大尺寸等离子体的粒子模拟方法，其特征在于，第四步中：

维持缩比前后模型中放电特征不改变，则必须要保证克努森数与霍尔参数在缩比前后模型中保持一致即均为一个常数：

式中，const表示常数，λ_ij为粒子碰撞自由程，ω_c为电子回旋频率，υ_ij为等离子体中i与j粒子之间碰撞频率；

进而得到：

式中，n为粒子数密度，σ为等离子体中i与j粒子之间的碰撞截面；

缩比前后模型中碰撞截面σ保持不变，缩比前后模型中粒子动能E_th以及粒子数密度n做相应变化，缩比前后模型中碰撞截面σ、粒子动能E_th以及粒子数密度n的缩比关系为：

式中，E_th为缩比模型中碰撞截面对应的粒子动能。

6.根据权利要求5所述的用于高密度、大尺寸等离子体的粒子模拟方法，其特征在于，第五步中：

缩比前后模型中推力器放电电流的变化关系为：

I_i，e～L^3/2

缩比前后模型中推力器推力和推进效率的变化关系为：

F～L²，η～1。