CN105678002A - 等离子体粒子-场自洽系统长期大规模高保真模拟方法 - Google Patents

Abstract

本发明公开了一种等离子体粒子-场自洽系统长期大规模高保真模拟方法，其包括：确定需要模拟计算的等离子体系统，获得相关表达方式；对所述相关表达方式进行离散处理，并根据离散处理后的表达方式离散Marsden-Weinstein泊松括号与系统的哈密顿量；根据离散后的Marsden-Weinstein泊松括号与系统的哈密顿量，获得等离子体粒子与电磁场的演化方程；利用欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，通过迭代求解获得离散格式，从而实现等离子体粒子-场自洽系统长期大规模高保真模拟。通过采用本发明公开的方法可以确保模拟过程的长期数值精确性与守恒性，并可确保高保真长期模拟。

Description

等离子体粒子-场自洽系统长期大规模高保真模拟方法

技术领域

本发明涉及等离子体数值模拟技术领域，尤其涉及一种等离子体粒子-场自洽系统长期大规模高保真模拟方法。

背景技术

等离子体科学与技术发展的目标是研究大量带电粒子与电磁场构成的自洽系统(即Vlasov-Maxwell，弗拉索夫-麦克斯韦系统)的演化与特性，并对其加以利用。等离子体是一种自由度极高的多尺度、强关联、具有集体效应的自洽系统，具有极高的复杂性。解析理论方法一般通过各种简化手段来探索等离子体的个别性质，以获得自洽系统的某个侧面信息。但解析方法很难处理普遍存在的多尺度、强关联、非线性等复杂过程。数值模拟方式通过将连续系统近似为离散系统，利用日益发展的强大计算资源，可以处理大量的、应用范围广泛的复杂等离子体问题。只要使用足够多的计算资源，通过数值模拟方法可以获得丰富的等离子体系统信息。因此，数值模拟方式在当前等离子体科学研究与工程应用中占据极其重要的地位。

在等离子体数值模拟领域，基于不同模型可设计各种不同模拟方法，其中PIC(Particle-in-Cell,质点网格法)方法作为一种基于第一性原理(弗拉索夫-麦克斯韦方程组)的数值模拟技术，可以还原等离子体系统的所有信息。因此PIC方法在等离子体科学研究与工业应用中已经得到广泛运用。

传统的PIC方法的实施思路是：

1、将电磁场、电流等场量按一定规则分配到空间网格上，例如利用计算电磁场常用Yee网格；同时将表示带电粒子的“宏粒子”采样点散布于空间网格中。

2、根据运动方程计算宏粒子在电磁场作用下在下一个时刻的位置和速度。运动方程一般为牛顿方程，所受的力由粒子位置处的电磁场计算，而某个位置的电磁场由格点上的场进行插值得到。

3、通过将宏粒子的运动坐标求和，计算格点上的电荷或电流分布。

4、利用格点上的电荷和电流分布，通过与麦克斯韦方程组对应的差分方程组计算下一个时刻位于格点上的电磁场值。

5、重复2-4步实现循环求解过程，可以得到最终时刻粒子-场自洽系统的状态和整个等离子体系统的完整演化过程。

从以上的模拟思路可以看出，传统PIC方法直接利用了格点上场与宏粒子的递推关系，并没有严格论证这种递推过程与原弗拉索夫-麦克斯韦系统的依赖关系和长期数值精度问题。最为重要的是，这种“粗糙”的离散递推方法不能保持系统的长期守恒性质(即模拟过程中的能量守恒、电荷守恒等)。实际计算中，数值误差的相干积累会不可避免地导致数值能量耗散和数值加热。由于数值误差的不断积累放大，在模拟多尺度、非线性等需要长期演化的等离子体过程时会导致结果严重失真，使模拟技术失效。可见，传统的PIC方法并不适用于解决多时间尺度、非线性物理问题，因而无法广泛应用于相应的等离子体科学研究与工程技术。

发明内容

本发明的目的是提供一种等离子体粒子-场自洽系统长期大规模高保真模拟方法，可以确保模拟过程的长期数值精确性与守恒性，并可确保高保真长期模拟。

本发明的目的是通过以下技术方案实现的：

一种等离子体粒子-场自洽系统长期大规模高保真模拟方法，包括：

确定需要模拟计算的等离子体系统，获得相关表达方式；

对所述相关表达方式进行离散处理，并根据离散处理后的表达方式离散Marsden-Weinstein泊松括号与系统的哈密顿量；

根据离散后的Marsden-Weinstein泊松括号与系统的哈密顿量，获得等离子体粒子与电磁场的演化方程；

利用欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，通过迭代求解获得离散格式，从而实现等离子体粒子-场自洽系统长期大规模高保真模拟。

进一步的，确定需要模拟计算的等离子体系统，获得相关表达方式包括：

确定需要模拟计算的等离子体系统的初始条件，获得粒子的分布函数和电磁场的解析表达方式。

进一步的，所述对所述相关表达方式进行离散处理包括：

利用Klimontovich表象离散粒子的分布函数，其表示为：

$f(x,p,t)={\mathrm{\Σ}}_{i=1}^N\mathrm{\δ}(x-X_i)\mathrm{\δ}(p-P_i);$

式中，f(x,p,t)为分布函数的离散量，x、p、t分别表示空间位置、正则动量、时间；X_i与P_i分别表示第i个采样点粒子的位置和正则动量；

将电磁场量(A,Y)离散到空间网格上，离散后表示为：(A_J,Y_J)；其中，A、Y分别表示磁矢势与电场的相反数，A_J与Y_J为序号为J的网格点上的值。

进一步的，所述根据离散处理后的表达方式离散Marsden-Weinstein泊松括号包括：

连续的Marsden-Weinstein泊松括号为：

$\{F,G\}=\∫f{\left\{\frac{\mathrm{\δ}F}{\mathrm{\δ}f},\frac{\mathrm{\δ}G}{\mathrm{\δ}f}\right\}}_{xp}dxdp+\∫f\left\{\frac{\mathrm{\δ}F}{\mathrm{\δ}A}\frac{\mathrm{\δ}G}{\mathrm{\δ}Y}-\frac{\mathrm{\δ}G}{\mathrm{\δ}A}\frac{\mathrm{\δ}F}{\mathrm{\δ}Y}\right\}dx;$

式中，F与G为任意物理量，f为离散后的分布函数；

离散后得到：

${\{F,G\}}_d=\underset{i=1}{\overset{N}{\Σ}}\left\{\frac{\mathrm{\δ}F}{{\mathrm{\δX}}_i}\frac{\mathrm{\δ}G}{{\mathrm{\δP}}_i}-\frac{\mathrm{\δ}G}{{\mathrm{\δX}}_i}\frac{\mathrm{\δ}F}{{\mathrm{\δP}}_i}\right\}+\underset{J=1}{\overset{M}{\Σ}}\left\{\frac{\mathrm{\δ}F}{{\mathrm{\δA}}_J}\frac{\mathrm{\δ}G}{{\mathrm{\δY}}_J}-\frac{\mathrm{\δ}G}{{\mathrm{\δA}}_J}\frac{\mathrm{\δ}F}{{\mathrm{\δY}}_J}\right\};$

式中，N为采样点粒子总数，M为网格点总数，{}_d表示离散泊松括号。

进一步的，离散系统的哈密顿量包括：

连续的系统的哈密度量为H(f,A,Y)：

$H(f,A,Y)=\frac{1}{2}\∫{(p-A)}^{2}fdxdp+\frac{1}{2}\∫\[Y^2+{(\▿\×A)}^2\]dx;$

其中,▽×表示旋度操作，上式右边第一项为粒子动能，第二项为场能；

离散后的系统的哈密顿量为：

$\begin{array}{c}\tilde{H}(X_i,P_i,A_J,Y_J)=\frac{1}{2}\underset{i=1}{\overset{N}{\Σ}}\[P_i^2-2P_i\·\underset{J=1}{\overset{M}{\Σ}}{A}_{J}W(X_i-x_J)\\ +\underset{J,L=1}{\overset{M}{\Σ}}{A}_J\·A_{L}W(X_i-x_J)W(X_i-x_L)\]+\frac{1}{2}\underset{J=1}{\overset{M}{\Σ}}\[Y_J^2+{({\▿}_d\×A_J)}^2\]\end{array};$

式中，L为网格点序号，W(X_i-x_J)表示任意一种差值函数，x_J、x_L分别为序号为J、L的网格点位置。

进一步的，根据离散后的Marsden-Weinstein泊松括号与系统的哈密顿量，获得等离子体粒子与电磁场的演化方程，该演化方程表示为：

${\stackrel{\·}{X}}_i={\{X_i,\tilde{H}\}}_d=P_i-\underset{J=1}{\overset{M}{\Σ}}{A}_{J}W(X_i-x_J);$

${\stackrel{\·}{A}}_J={\{A_J,\tilde{H}\}}_d=Y_J;$

${\stackrel{\·}{P}}_i={\{P_i,\tilde{H}\}}_d=\underset{J=1}{\overset{M}{\Σ}}(P_i\·A_J)\▿W(X_i-x_J)-\underset{J,L=1}{\overset{M}{\Σ}}(A_J\·A_L)W(X_i-x_J)\▿W(X_i-x_L);$

${\stackrel{\·}{Y}}_J={\{Y_J,\tilde{H}\}}_d=\underset{i=1}{\overset{N}{\Σ}}{P}_{i}W(X_i-x_J)-\underset{i=1}{\overset{N}{\Σ}}\underset{J,L=1}{\overset{M}{\Σ}}{A}_{L}W(X_i-x_J)W(X_i-x_L)-{({\▿}_d\×{\▿}_d\×A)}_J;$

式中，▽_d×为离散旋度符号，则有：

${({\▿}_d\×A)}_J=\left(\begin{array}{c}\frac{A_{i^{\′},j,k}^3-A_{i^{\′},j-1,k}^3}{\mathrm{\Δ}y}-\frac{A_{i^{\′},j,k}^2-A_{i^{\′},j,k-1}^2}{\mathrm{\Δ}z}\\ \frac{A_{i^{\′},j,k}^1-A_{i^{\′},j,k-1}^1}{\mathrm{\Δ}z}-\frac{A_{i^{\′},j,k}^3-A_{i^{\′}-1,j,k}^3}{\mathrm{\Δ}x}\\ \frac{A_{i^{\′},j,k}^2-A_{i^{\′}-1,j,k}^2}{\mathrm{\Δ}z}-\frac{A_{i^{\′},j,k}^1-A_{i^{\′},j-1,k}^1}{\mathrm{\Δ}x}\end{array}\right);$

其中，J＝(i',j,k)，Δx、Δy、Δz分别表示x、y、z方向的网格长度；为离散磁矢势的三个分量。

进一步的，所述利用欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，通过迭代求解获得离散格式，从而实现等离子体粒子-场自洽系统长期大规模高保真模拟包括：

欧拉-辛算法表示为：

$p^{n+1}=p^n-\mathrm{\Δ}t\frac{\∂H}{\∂x}(p^{n+1},x^n);$

$x^{n+1}=x^n-\mathrm{\Δ}t\frac{\∂H}{\∂p}(p^{n+1},x^n);$

式中，n表示计算的步数，Δt表示计算时的时间步长；；

利用所述欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，表示为：

$\frac{X_i^{n+1}-X_i^n}{\mathrm{\Δ}t}=P_i^{n+1}-\underset{J=1}{\overset{M}{\Σ}}{A}_J^{n}W(X_i^n-x_J);$

$\frac{A_J^{n+1}-A_J^n}{\mathrm{\Δ}t}=Y_J^{n+1};$

$\frac{P_i^{n+1}-P_i^n}{\mathrm{\Δ}t}=\underset{J=1}{\overset{M}{\Σ}}(P_i^{n+1}\·A_J^n)\▿W(X_i^n-x_J)-\underset{J,L=1}{\overset{M}{\Σ}}(A_J^n\·A_L^n)W(X_i^n-x_J)\▿W(X_i^n-x_L);$

$\frac{Y_J^{n+1}-Y_J^n}{\mathrm{\Δ}t}=\underset{i=1}{\overset{N}{\Σ}}{P}_i^{n+1}W(X_i^n-x_J)-\underset{i=1}{\overset{N}{\Σ}}\underset{J,L=1}{\overset{M}{\Σ}}{A}_L^{n}W(X_i^n-x_J)W(X_i^n-x_L)-{({\▿}_d\×{\▿}_d\×A^n)}_J;$

通过设定电磁场与粒子的初始条件与边界条件，并利用上述表达方式进行时间步进，最终实现等离子体粒子-场自洽系统长期大规模高保真模拟。

由上述本发明提供的技术方案可以看出，首先，通过特殊离散几何微分格式离散了空间电磁场，保证了电磁场模拟的自洽性与守恒律；其次，通过离散泊松括号和离散哈密顿量构造离散哈密顿系统，实现辛PIC方法，保证了模拟过程的长期数值精确性与守恒性，从而保证了高保真长期模拟；同时，在对带电粒子分布离散过程中应用了光滑因子以消除数值噪音的影响，解决了PIC模拟中的数值噪音问题，使物理信号更容易识别；并且，利用半隐式递推关系使得每次步进解的过程局域化，有利于并行化计算的实现；还通过进程-线程混合并行技术，实现高保真正则辛PIC方法长期大规模并行模拟的高并行扩展性；另外，由于本方案能够长期保证模拟的精度和守恒性，可应用于多尺度、非线性、强关联等各种情况，并具有极高的大规模计算效率，能够灵活、高效地广泛应用于各种等离子体系统和过程，具有极高的应用和扩展价值。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的一种等离子体粒子-场自洽系统长期大规模高保真模拟方法的流程图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

图1为本发明实施例提供的一种等离子体粒子-场自洽系统长期大规模高保真模拟方法的流程图。如图1所示，其主要包括如下步骤：

步骤11、确定需要模拟计算的等离子体系统，获得相关表达方式。

本发明实施例中，针对需要模拟计算的等离子体系统，可以采用常规方法来确定需要模拟计算的等离子体系统的初始条件，从而获得粒子的分布函数和电磁场的解析表达方式。

步骤12、对所述相关表达方式进行离散处理，并根据离散处理后的表达方式离散Marsden-Weinstein泊松括号与系统的哈密顿量。

本步骤主要包括三个处理过程：

一、离散化分布函数和电磁场。

本发明实施例中，利用Klimontovich表象离散粒子的分布函数，其表示为：

$f(x,p,t)={\mathrm{\Σ}}_{i=1}^N\mathrm{\δ}(x-X_i)\mathrm{\δ}(p-P_i);$

式中，f(x,p,t)为分布函数的离散量，x、p、t分别表示空间位置、正则动量、时间；X_i与P_i分别表示第i个采样点粒子的位置和正则动量；

同时，将电磁场量(A,Y)离散到空间网格上，离散后表示为：(A_J,Y_J)；其中，A、Y分别表示磁矢势与电场的相反数，A_J与Y_J为序号为J的网格点上的值。

通过采用Klimontovich表象离散分布函数，将连续的多元函数转化为有限个正则采样，极大地减少了计算难度；在对带电粒子分布离散过程中应用了光滑因子以消除数值噪音的影响，解决了PIC模拟中的数值噪音问题，使物理信号更容易识别。

二、离散Marsden-Weinstein泊松括号

本发明实施例中，连续的Marsden-Weinstein泊松括号为：

$\{F,G\}=\∫f{\left\{\frac{\mathrm{\δ}F}{\mathrm{\δ}f},\frac{\mathrm{\δ}G}{\mathrm{\δ}f}\right\}}_{xp}dxdp+\∫f\left\{\frac{\mathrm{\δ}F}{\mathrm{\δ}A}\frac{\mathrm{\δ}G}{\mathrm{\δ}Y}-\frac{\mathrm{\δ}G}{\mathrm{\δ}A}\frac{\mathrm{\δ}F}{\mathrm{\δ}Y}\right\}dx;$

式中，F与G为任意物理量，f为离散后的分布函数；

利用前文离散后的分布函数和电磁场离散上述Marsden-Weinstein泊松括号，得到：

${\{F,G\}}_d=\underset{i=1}{\overset{N}{\Σ}}\left\{\frac{\mathrm{\δ}F}{{\mathrm{\δX}}_i}\frac{\mathrm{\δ}G}{{\mathrm{\δP}}_i}-\frac{\mathrm{\δ}G}{{\mathrm{\δX}}_i}\frac{\mathrm{\δ}F}{{\mathrm{\δP}}_i}\right\}+\underset{J=1}{\overset{M}{\Σ}}\left\{\frac{\mathrm{\δ}F}{{\mathrm{\δA}}_J}\frac{\mathrm{\δ}G}{{\mathrm{\δY}}_J}-\frac{\mathrm{\δ}G}{{\mathrm{\δA}}_J}\frac{\mathrm{\δ}F}{{\mathrm{\δY}}_J}\right\};$

式中，N为采样点粒子总数，M为网格点总数，{}_d表示离散泊松括号。

本发明实施例中，采用连续的Marsden-Weinstein泊松括号体现了物理系统结构的本质；通过直接离散泊松括号，保证不破坏原连续系统的物理结构，达成对原系统的最好的离散描述。

三、离散系统的哈密顿量

所述系统即为自洽系统(Vlasov-Maxwell系统)，已知连续的Vlasov-Maxwell系统的哈密顿量为H(f,A,Y)：

$H(f,A,Y)=\frac{1}{2}\∫{(p-A)}^{2}fdxdp+\frac{1}{2}\∫\[Y^2+{(\▿\×A)}^2\]dx;$

其中,▽×表示旋度操作,上式右边第一项为粒子动能，第二项为场能；

通过插值离散场点的值得到粒子位置的电磁场，即：

$A(X_i,t)=\underset{J=1}{\overset{M}{\Σ}}{A}_{J}W(X_i-x_J)$

$Y(X_i,t)=\underset{J=1}{\overset{M}{\Σ}}{Y}_J\left(t\right)W(X_i-x_J);$

则离散后的系统的哈密顿量为：

$\begin{array}{c}\tilde{H}(X_i,P_i,A_J,Y_J)=\frac{1}{2}\underset{i=1}{\overset{N}{\Σ}}\[P_i^2-2P_i\·\underset{J=1}{\overset{M}{\Σ}}{A}_{J}W(X_i-x_J)\\ +\underset{J,L=1}{\overset{M}{\Σ}}{A}_J\·A_{L}W(X_i-x_J)W(X_i-x_L)\]+\frac{1}{2}\underset{J=1}{\overset{M}{\Σ}}\[Y_J^2+{({\▿}_d\×A_J)}^2\]\end{array};$

式中，L为网格点序号，W(X_i-x_J)表示任意一种差值函数，x_J、x_L分别为序号为J、L的网格点位置。

本发明实施例中，通过对正则哈密顿量的离散实现了系统正则辛结构的保持，为长期数值稳定性与精确性和辛欧拉方法的实施提供了保障。

步骤13、根据离散后的Marsden-Weinstein泊松括号与系统的哈密顿量，获得等离子体粒子与电磁场的演化方程。

利用前述步骤12离散获得的泊松括号与系统的哈密顿量得到每个粒子和电磁场的演化方程，即任意一个物理量F的演化可以由给出。

则该演化方程可以表示为：

${\stackrel{\·}{X}}_i={\{X_i,\tilde{H}\}}_d=P_i-\underset{J=1}{\overset{M}{\Σ}}{A}_{J}W(X_i-x_J);$

${\stackrel{\·}{A}}_J={\{A_J,\tilde{H}\}}_d=Y_J;$

${\stackrel{\·}{P}}_i={\{P_i,\tilde{H}\}}_d=\underset{J=1}{\overset{M}{\Σ}}(P_i\·A_J)\▿W(X_i-x_J)-\underset{J,L=1}{\overset{M}{\Σ}}(A_J\·A_L)W(X_i-x_J)\▿W(X_i-x_L);$

${\stackrel{\·}{Y}}_J={\{Y_J,\tilde{H}\}}_d=\underset{i=1}{\overset{N}{\Σ}}{P}_{i}W(X_i-x_J)-\underset{i=1}{\overset{N}{\Σ}}\underset{J,L=1}{\overset{M}{\Σ}}{A}_{L}W(X_i-x_J)W(X_i-x_L)-{({\▿}_d\×{\▿}_d\×A)}_J;$

式中，▽_d×为离散旋度符号，则有：

${({\▿}_d\×A)}_J=\left(\begin{array}{c}\frac{A_{i^{\′},j,k}^3-A_{i^{\′},j-1,k}^3}{\mathrm{\Δ}y}-\frac{A_{i^{\′},j,k}^2-A_{i^{\′},j,k-1}^2}{\mathrm{\Δ}z}\\ \frac{A_{i^{\′},j,k}^1-A_{i^{\′},j,k-1}^1}{\mathrm{\Δ}z}-\frac{A_{i^{\′},j,k}^3-A_{i^{\′}-1,j,k}^3}{\mathrm{\Δ}x}\\ \frac{A_{i^{\′},j,k}^2-A_{i^{\′}-1,j,k}^2}{\mathrm{\Δ}z}-\frac{A_{i^{\′},j,k}^1-A_{i^{\′},j-1,k}^1}{\mathrm{\Δ}x}\end{array}\right);$

其中，J＝(i',j,k)，本发明实施例中，(i',j,k)表示(x,y,z)三个方向的序号；Δx、Δy、Δz分别表示x、y、z方向的网格长度；为离散磁矢势的三个分量。

本发明实施例中，将泊松括号作用于正则哈密顿量构造系统的演化方程；在保证原系统守恒量的前提下，极大程度地减少计算量。

步骤14、利用欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，通过迭代求解获得离散格式，从而实现等离子体粒子-场自洽系统长期大规模高保真模拟。

本发明实施例中，使用欧拉-辛算法处理步骤13得到系统的演化方程，得到正则辛PIC算法，使用得到的正则辛PIC算法进行等离子体模拟计算。在本实施例中，欧拉-辛算法表示为：

$p^{n+1}=p^n-\mathrm{\Δ}t\frac{\∂H}{\∂x}(p^{n+1},x^n);$

$x^{n+1}=x^n-\mathrm{\Δ}t\frac{\∂H}{\∂p}(p^{n+1},x^n);$

式中，n表示计算的步数，例如，pⁿ⁺¹表示第n+1步计算时的正则动量，xⁿ⁺¹表示第n+1步计算时的空间位置；Δt表示计算时的时间步长；

利用所述欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，表示为：

$\frac{X_i^{n+1}-X_i^n}{\mathrm{\Δ}t}=P_i^{n+1}-\underset{J=1}{\overset{M}{\Σ}}{A}_J^{n}W(X_i^n-x_J);$

$\frac{A_J^{n+1}-A_J^n}{\mathrm{\Δ}t}=Y_J^{n+1};$

$\frac{P_i^{n+1}-P_i^n}{\mathrm{\Δ}t}=\underset{J=1}{\overset{M}{\Σ}}(P_i^{n+1}\·A_J^n)\▿W(X_i^n-x_J)-\underset{J,L=1}{\overset{M}{\Σ}}(A_J^n\·A_L^n)W(X_i^n-x_J)\▿W(X_i^n-x_L);$

$\frac{Y_J^{n+1}-Y_J^n}{\mathrm{\Δ}t}=\underset{i=1}{\overset{N}{\Σ}}{P}_i^{n+1}W(X_i^n-x_J)-\underset{i=1}{\overset{N}{\Σ}}\underset{J,L=1}{\overset{M}{\Σ}}{A}_L^{n}W(X_i^n-x_J)W(X_i^n-x_L)-{({\▿}_d\×{\▿}_d\×A^n)}_J;$

本发明实施例中，通过辛欧拉保证时间步进格式的辛结构，从而保证模拟的长期守恒性质与高保真性质；同时，利用半隐式递推关系使得每次步进解的过程局域化，有利于并行化计算的实现。

最后，通过设定电磁场与粒子的初始条件与边界条件，并利用上述表达方式进行时间步进，最终实现等离子体粒子-场自洽系统长期大规模高保真模拟。

本发明实施例的上述方案相对于传统方案而言，主要具有如下优点：

1)通过特殊离散几何微分格式离散了空间电磁场，保证了电磁场模拟的自洽性与守恒律。

2)通过离散泊松括号和离散哈密顿量构造离散哈密顿系统，实现辛PIC方法，保证了模拟过程的长期数值精确性与守恒性，从而保证了高保真长期模拟。

3)在对带电粒子分布离散过程中应用了光滑因子以消除数值噪音的影响，解决了PIC模拟中的数值噪音问题，使物理信号更容易识别。

4)利用半隐式递推关系使得每次步进解的过程局域化，有利于并行化计算的实现。

5)通过进程-线程混合并行技术，实现高保真正则辛PIC方法长期大规模并行模拟的高并行扩展性。

6)由于本方法能够长期保证模拟的精度和守恒性，可应用于多尺度、非线性、强关联等各种情况，并具有极高的大规模计算效率，能够灵活、高效地广泛应用于各种等离子体系统和过程，具有极高的应用和扩展价值。

为了进一步说明本发明，下面结合一具体示例做详细说明。需要强调的是，下述示例中所涉及的具体数值仅为举例，并非构成限制。在实际应用中，各个参数的具体数值可根据实际情况来设定。

本示例中，具体处理步骤如下：

第一步：确定需要模拟计算的等离子体系统；

针对朗道阻尼过程的问题，朗道阻尼问题采用一维模型描述，关于连续分布函数，粒子在空间是均匀分布，速度空间满足麦克斯韦分布，没有外加电磁场，只是在系统空间加一个扰动电场，以便观察其阻尼情况，本问题是一个严格的Vlasov-Maxwell系统的问题。采用现有技术确定等离子体物理系统的初始条件，包括粒子的分布函数和电磁场的解析表达方式。

第二步：离散分布函数和电磁场；

利用Klimontovich表象离散第一步得到的分布函数，同时将电磁场(A,Y)离散到空间网格上。离散过程中，将空间划分为960×2×2的方网格，网格的长度设为Δx＝2.4535×10^-4m，所有的场值都只在格点上取值，即离散场；求解的电子密度是n_e＝1.239×10¹⁹/m³，每个网格分布10³个电子采样点，按照热速度为v_T＝0.2c的麦克斯韦分布散入粒子，通过采样点即实现了分布函数的离散，这样采样粒子的总数为N＝9.57×10⁷，自由度数目为D＝5.74×10⁸，这实现了分布函数的离散。

第三步：离散Marsden-Weinstein泊松括号；

通过第一步离散得到的分布函数和电磁场，将连续的Marsden-Weinstein泊松括号离散。具体公式可参见前述步骤12，本示例中仅离散x方向的分量。

第四步：离散Vlasov-Maxwell系统的哈密顿量；

离散Vlasov-Maxwell系统的哈密顿量。将离散哈密顿量的所有矢量换为其x方向的分量，另外，格点值相应的简化为一维格点，最终得到的离散格式关于格点的项都可以相应的简化，即J的值变为0-3839的序号。最后，插值函数选择6阶差值函数，即：

$W\left(q\right)=\left\{\begin{array}{cc}0,& q>2\\ \frac{30}{128}{q}^6-\frac{27}{32}{q}^5+\frac{33}{64}{q}^4-q+1,& 1<q\&le;2\\ \frac{7}{16}{q}^6-\frac{23}{32}{q}^5+\frac{175}{256}{q}^4-\frac{109}{128}{q}^2+\frac{399}{512},& 0<q\&le;1\\ \frac{9}{16}{q}^6+\frac{19}{32}{q}^5+\frac{175}{256}{q}^4-\frac{109}{128}{q}^2+\frac{399}{512},& -1<q\&le;0\\ \frac{47}{128}{q}^6+\frac{23}{32}{q}^5+\frac{37}{64}{q}^4+q+1,& -2<q\&le;-1\\ 0,& q<-2\end{array}\right.$

其中，q＝X_i-x_J；

第五步：利用离散哈密顿量和泊松括号获得粒子、场的演化方程；

利用第三步离散得到的泊松括号和第四步离散得到的哈密顿量得到每个粒子和离散场的演化方程。在本实施例中，获得的离散系统的演化方程同样只包含x分量，原来的3维适量方程相应的简化。

第六步：利用欧拉-辛算法处理得到离散系统的演化方程；

采用欧拉-辛算法形式处理第五步得到离散系统方程，即得离散格式。具体过程可参见前文的步骤14。

第七步：迭代求解第五步得到的离散格式；

通过设定场与粒子的初始条件与边界条件，利用第二步到第六步得到的高保真正则辛PIC格式进行时间步进，实现等离子体系统第一性原理的大规模高保真模拟。计算时，时间步长为Δt＝Δx/2c。初始电场扰动为E₁＝E₁cos(kx)e_x，其中波矢为k＝2π/224Δx，扰动场振幅为E₁＝9.103×10⁴v/m。通过模拟80000时间步，可以揭示非线性朗道阻尼的完整图像。正如使用辛算法所期望的，能量的数值误差不会增长并且始终局域在百万分之一的范围内。利用色散关系得到的理论阻尼率是ω_i＝-1.3926×10⁹/s，理论实频率是ω_r＝9.116×10⁹/s。

经过验证，应用本发明方法的模拟与理论符合度很高，经过10⁷步数的计算，模拟结果与理论计算的最大误差不超过10^-6，且不随计算步数增加而累积，与现有方法的误差随计算步数增加累积相比，优势极为明显。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (7)

1.一种等离子体粒子-场自洽系统长期大规模高保真模拟方法，其特征在于，包括：

确定需要模拟计算的等离子体系统，获得相关表达方式；

对所述相关表达方式进行离散处理，并根据离散处理后的表达方式离散Marsden-Weinstein泊松括号与系统的哈密顿量；

根据离散后的Marsden-Weinstein泊松括号与系统的哈密顿量，获得等离子体粒子与电磁场的演化方程；

利用欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，通过迭代求解获得离散格式，从而实现等离子体粒子-场自洽系统长期大规模高保真模拟。

2.根据权利要求1所述的方法，其特征在于，确定需要模拟计算的等离子体系统，获得相关表达方式包括：

确定需要模拟计算的等离子体系统的初始条件，获得粒子的分布函数和电磁场的解析表达方式。

3.根据权利要求1或2所述的方法，其特征在于，所述对所述相关表达方式进行离散处理包括：

利用Klimontovich表象离散粒子的分布函数，其表示为：

$f(x,p,t)={\mathrm{\&Sigma;}}_{i=1}^N\mathrm{\&delta;}(x-X_i)\mathrm{\&delta;}(p-P_i);$

式中，f(x,p,t)为分布函数的离散量，x、p、t分别表示空间位置、正则动量、时间；X_i与P_i分别表示第i个采样点粒子的位置和正则动量；

将电磁场量(A,Y)离散到空间网格上，离散后表示为：(A_J,Y_J)；其中，A、Y分别表示磁矢势与电场的相反数，A_J与Y_J为序号为J的网格点上的值。

4.根据权利要求3所述的方法，其特征在于，所述根据离散处理后的表达方式离散Marsden-Weinstein泊松括号包括：

连续的Marsden-Weinstein泊松括号为：

$\{F,G\}=\&Integral;f{\{\frac{\mathrm{\&delta;}F}{\mathrm{\&delta;}f},\frac{\mathrm{\&delta;}G}{\mathrm{\&delta;}f}\}}_{xp}dxdp+\&Integral;f\{\frac{\mathrm{\&delta;}F}{\mathrm{\&delta;}A}\frac{\mathrm{\&delta;}G}{\mathrm{\&delta;}Y}-\frac{\mathrm{\&delta;}G}{\mathrm{\&delta;}A}\frac{\mathrm{\&delta;}F}{\mathrm{\&delta;}Y}\}dx;$

式中，F与G为任意物理量，f为离散后的分布函数；

离散后得到：

${\{F,G\}}_d=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\{\frac{\mathrm{\&delta;}F}{{\mathrm{\&delta;X}}_i}\frac{\mathrm{\&delta;}G}{{\mathrm{\&delta;P}}_i}-\frac{\mathrm{\&delta;}G}{{\mathrm{\&delta;X}}_i}\frac{\mathrm{\&delta;}F}{{\mathrm{\&delta;P}}_i}\}+\underset{J=1}{\overset{M}{\mathrm{\&Sigma;}}}\{\frac{\mathrm{\&delta;}F}{{\mathrm{\&delta;A}}_J}\frac{\mathrm{\&delta;}G}{{\mathrm{\&delta;Y}}_J}-\frac{\mathrm{\&delta;}G}{{\mathrm{\&delta;A}}_J}\frac{\mathrm{\&delta;}F}{{\mathrm{\&delta;Y}}_J}\};$

式中，N为采样点粒子总数，M为网格点总数，{}_d表示离散泊松括号。

5.根据权利要求4所述的方法，其特征在于，离散系统的哈密顿量包括：

连续的系统的哈密度量为H(f,A,Y)：

$H(f,A,Y)=\frac{1}{2}\&Integral;{(p-A)}^{2}fdxdp+\frac{1}{2}\&Integral;\&lsqb;Y^2+{(\&dtri;\&times;A)}^2\&rsqb;dx;$

其中,表示旋度操作，上式右边第一项为粒子动能，第二项为场能；

离散后的系统的哈密顿量为：

$\begin{array}{c}\tilde{H}(X_i,P_i,A_J,Y_J)=\frac{1}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\&lsqb;P_i^2-2P_i\&CenterDot;\underset{J=1}{\overset{M}{\mathrm{\&Sigma;}}}{A}_{J}W(X_i-x_J)\\ +\underset{J,L=1}{\overset{M}{\mathrm{\&Sigma;}}}{A}_J\&CenterDot;A_{L}W(X_i-x_J)W(X_i-x_L)\&rsqb;+\frac{1}{2}\underset{J=1}{\overset{M}{\mathrm{\&Sigma;}}}\&lsqb;Y_J^2+{({\&dtri;}_d\&times;A_J)}^2\&rsqb;\end{array};$

式中，L为网格点序号，W(X_i-x_J)表示任意一种差值函数，x_J、x_L分别为序号为J、L的网格点位置。

6.根据权利要求5所述的方法，其特征在于，根据离散后的Marsden-Weinstein泊松括号与系统的哈密顿量，获得等离子体粒子与电磁场的演化方程，该演化方程表示为：

${\stackrel{\&CenterDot;}{X}}_i={\{X_i,\tilde{H}\}}_d=P_i-\underset{J=1}{\overset{M}{\&Sigma;}}{A}_{J}W(X_i-x_J);$

${\stackrel{\&CenterDot;}{A}}_J={\{A_J,\tilde{H}\}}_d=Y_J;$

${\stackrel{\&CenterDot;}{P}}_i={\{P_i,\tilde{H}\}}_d=\underset{J=1}{\overset{M}{\&Sigma;}}(P_i\&CenterDot;A_J)\&dtri;W(X_i-x_J)-\underset{J,L=1}{\overset{M}{\&Sigma;}}(A_J\&CenterDot;A_L)W(X_i-x_J)\&dtri;W(X_i-x_L);$

${\stackrel{\&CenterDot;}{Y}}_J={\{Y_J,\tilde{H}\}}_d=\underset{i=1}{\overset{N}{\&Sigma;}}{P}_{i}W(X_i-x_J)-\underset{i=1}{\overset{N}{\&Sigma;}}\underset{J,L=1}{\overset{M}{\&Sigma;}}{A}_{L}W(X_i-x_J)W(X_i-x_L)-{({\&dtri;}_d\&times;{\&dtri;}_d\&times;A)}_J;$

式中，为离散旋度符号，则有：

${({\&dtri;}_d\&times;A)}_J=\left(\begin{array}{c}\frac{A_{i^{\&prime;},j,k}^3-A_{i^{\&prime;},j-1,k}^3}{\mathrm{\&Delta;}y}-\frac{A_{i^{\&prime;},j,k}^2-A_{i^{\&prime;},j,k-1}^2}{\mathrm{\&Delta;}z}\\ \frac{A_{i^{\&prime;},j,k}^1-A_{i^{\&prime;},j,k-1}^1}{\mathrm{\&Delta;}z}-\frac{A_{i^{\&prime;},j,k}^3-A_{i^{\&prime;}-1,j,k}^3}{\mathrm{\&Delta;}x}\\ \frac{A_{i^{\&prime;},j,k}^2-A_{i^{\&prime;}-1,j,k}^2}{\mathrm{\&Delta;}z}-\frac{A_{i^{\&prime;},j,k}^1-A_{i^{\&prime;},j-1,k}^1}{\mathrm{\&Delta;}x}\end{array}\right);$

其中，J＝(i',j,k)，Δx、Δy、Δz分别表示x、y、z方向的网格长度；为离散磁矢势的三个分量。

7.根据权利要求6所述的方法，其特征在于，所述利用欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，通过迭代求解获得离散格式，从而实现等离子体粒子-场自洽系统长期大规模高保真模拟包括：

欧拉-辛算法表示为：

$p^{n+1}=p^n-\mathrm{\&Delta;}t\frac{\&part;H}{\&part;x}(p^{n+1},x^n);$

$x^{n+1}=x^n-\mathrm{\&Delta;}t\frac{\&part;H}{\&part;p}(p^{n+1},x^n);$

式中，n表示计算的步数，Δt表示计算时的时间步长；；

利用所述欧拉-辛算法处理所述等离子体粒子与电磁场的演化方程，表示为：

$\frac{X_i^{n+1}-X_i^n}{\mathrm{\&Delta;}t}=P_i^{n+1}-\underset{J=1}{\overset{M}{\&Sigma;}}{A}_J^{n}W(X_i^n-x_J);$

$\frac{A_J^{n+1}-A_J^n}{\mathrm{\&Delta;}t}=Y_J^{n+1};$

$\frac{P_i^{n+1}-P_i^n}{\mathrm{\&Delta;}t}=\underset{J=1}{\overset{M}{\&Sigma;}}(P_i^{n+1}\&CenterDot;A_J^n)\&dtri;W(X_i^n-x_J)-\underset{J,L=1}{\overset{M}{\&Sigma;}}(A_J^n\&CenterDot;A_L^n)W(X_i^n-x_J)\&dtri;W(X_i^n-x_L);$

$\frac{Y_J^{n+1}-Y_J^n}{\mathrm{\&Delta;}t}=\underset{i=1}{\overset{N}{\&Sigma;}}{P}_i^{n+1}W(X_i^n-x_J)-\underset{i=1}{\overset{N}{\&Sigma;}}\underset{J,L=1}{\overset{M}{\&Sigma;}}{A}_L^{n}W(X_i^n-x_J)W(X_i^n-x_L)-{({\&dtri;}_d\&times;{\&dtri;}_d\&times;A^n)}_J;$

通过设定电磁场与粒子的初始条件与边界条件，并利用上述表达方式进行时间步进，最终实现等离子体粒子-场自洽系统长期大规模高保真模拟。