电推力器中复杂电场求解的超松弛并行算法

，，*，，，，

1.大连理工大学 工业装备结构分析国家重点实验室，大连 116024 2.中国空间技术研究院 通信卫星事业部，北京 100094

电推力器中复杂电场求解的超松弛并行算法

韩亚杰1，夏广庆1，*，吴秋云1，陈留伟1，周念东1，温正2

1.大连理工大学 工业装备结构分析国家重点实验室，大连 116024 2.中国空间技术研究院 通信卫星事业部，北京 100094

目前，电推力器中复杂电场的三维泊松方程求解仍存在速度偏慢、并行性欠佳等缺陷。通过分解三维求解区域，针对复杂电场提出一种超松弛并行(P-SOR)迭代算法。将求解区域划分为多个子块，利用超松弛迭代格式构造出若干分组显式格式，分别给出不同迭代步数下的求解方程以进行并行计算。通过数值模拟，计算时间至少缩短一半以上，该P-SOR方法较传统迭代格式在快速性方面取得较大进展，对电推进领域的数值仿真研究起到促进作用。

电推力器；泊松方程；超松弛迭代；并行算法

泊松方程是等离子体物理、静电学中常见的偏微分方程，其应用范围较广，特别是在电推力器中复杂电场的求解方面。现阶段，大规模的科学工程计算对于速度快、容量大的并行机和有效的数值并行算法均有较大的需求，特别是在航天领域。基于粒子云网格(Particle In Cell，PIC)方法[1]的数值计算，通过泊松方程来求解复杂电场，在电帆仿真计算[2]，霍尔推力器的相关研究[3]等领域均有广泛应用。但由于三维模型求解速度很慢，诸多研究人员针对模型做出简化，如采用二维模型、二维轴对称模型[4-5]等，其精确度难免会下降。所以本文针对加快三维泊松方程的求解速度开展研究。

国内外已有一些专家学者在该领域进行了探索。现在通用的并行模块大多采用Jacobi迭代[6]，其具有显性的并行性。但是对于求解速度，超松弛迭代等方法则要远远高于Jacobi格式，也有很多学者对其进行研究。例如，文献[7]提出一种针对点的SOR并行方法，但是在扩展性方面存在较大的劣势，不利于大规模应用。文献[8]提出一种新型的、应用于统一计算架构(CUDA)的并行算法，利用GPU的硬件特性对泊松方程的求解进行加速。文献[9]提出一种块SOR理论，改进其求解速度，并对收敛性做出研究。文献[10]提出一种基于Richardson外推法的三维泊松方程差分格式，具有使用网格基架点少、精度高、稳定性好等优点。但是对于三维泊松方程的多线程并行计算的研究仍较为欠缺。

本文提出一种基于三维泊松方程的并行计算理论，适用于OpenMP并行模块[11],并采用电场结构较为复杂的柱形和球形惯性静电约束电推力器[12](Inertial Electrostatic Confinement Thruster，IECT)为模型进行相应的数值模拟，通过对比不同迭代方法的求解时间，对其快速性进行验证。

1 泊松问题及分组显式格式

1.1泊松问题

考虑在长方体区域Ω上的三维泊松方程，

式中：Ω=[0≤x≤L]×[0≤y≤M]×[0≤z≤N]为求解区域，∂Ω为Ω的边界；u(x,y,z)为待求函数，f(x,y,z)和g(x,y,z)为已知光滑函数。

对于泊松方程内部点的逼近方法有很多，这里采用最为常见的七点差分格式[13]：

对区域Ω进行一致剖分，空间步长取为h，在x，y和z方向分别为(l+1)h=L，(m+1)h=M与(n+1)h=N，则定义xi=ih,i=0,1,…,l+1；yj=jh,j=0,1,…,m+1；zk=kh,k=0,1,…,n+1，记(xi,yj,zk)为(i,j,k)，则有：

SOR方法很符合求解此类问题，但它没有明显的并行性，所以这里将给出一个新的并行SOR迭代算法，简称为P-SOR。根据对二维并行迭代格式的研究[14]，推导出基于三维泊松方程的8个基本SOR迭代格式：

式中：ω为松弛因子。

1.2 分组显式格式

首先，考虑求解域中相邻8个网格点的数值解，如图1所示，联立A1～A8可形成(8×8)阶方程组，采用分组显式的思想，即可得到求解域节点(i,j,k)，(i,j+1,k)，(i+1,j,k)，(i+1,j+1,k)，(i,j,k+1)，(i,j+1,k+1)，(i+1,j,k+1)，(i+1,j+1,k+1)处的数值解公式。

图1 中间点迭代示意Fig.1 Schematic diagram of intermediate point iteration

令

其中

即可求出方程组：

f2d5+4ω2d6+4ω2d7+2ω3d8)

4ω2d5+f2d6+2ω3d7+4ω2d8)

4ω2d5+2ω3d6+f2d7+4ω2d8)

2ω3d5+4ω2d6+4ω2d7+f2d8)

f1d5+f2d6+f2d7+4ω2d8)

f2d5+f1d6+4ω2d7+f2d8)

f2d5+4ω2d6+f1d7+f2d8)

4ω2d5+f2d6+f2d7+f1d8)

式中：

f1=-14ω2+72，f2=-ω3+12ω

类似地，利用分组显式基于“四点一组”的思想可构造其他6个新的格式。

利用格式A1,A3,A5,A7求解内边界点(i,j0,k)，(i+1,j0,k)，(i,j0,k+1)和(i+1,j0,k+1)，j0=j+2,j+3,…,m可得格式A10：

用类似的方法即可求出另外5个迭代格式A11～A15，其基本表现形式相同，这里不再详述。

特别地，当f(i,j,k)=0时，该泊松方程即退化为拉普拉斯方程。对于超松弛系数ω的选取，可以参考文献[15]。

2 P-SOR算法流程

2.1区域分解

首先，将区域Ω分解为若干子域，这里以8个非重叠子域为例(其他以此类推)，如图2所示，用6个网格平面(i=p,p+1,j=q,q+1,k=r,r+1)将求解区域划分的8个子域，其中p,q,r为任意正整数，0

图2 求解域网格划分Fig.2 Meshing partitions in the calculation domain

2.2 迭代算法

显然，在求解过程中对于相同的区域不能使用同一种迭代格式，这样会造成不同时间步的电势值无法更新迭代。所以，对于奇数次迭代和偶数次迭代，采用不同的迭代算法进行实现，步骤如下：

1)迭代次数为奇数。

利用格式A1～A8依次对子域Ω1～Ω8中的网格点由外向内进行迭代。

2)迭代次数为偶数。

利用格式A9求解网格中心交叉点(p,q,r)，(p+1,q,r)，(p,q+1,r)，(p+1,q+1,r)，(p,q,r+1)，(p+1,q,r+1)，(p,q+1,r+1)，(p+1,q+1,r+1)。

利用格式A10求解内边界点(p,j,r)，(p+1,j,r)，(p,j,r+1)和(p+1,j,r+1)，j=q+2,q+3,…,m；利用格式A11求解内边界点(p,j,r)，(p+1,j,r)，(p,j,r+1)和(p+1,j,r+1)，j=q-1,q-2,…,1；利用格式A12求解内边界点(i,q,r)，(i,q+1,r)，(i,q,r+1)和(i,q+1,r+1)，i=p-1,p-2,…,1；利用格式A13求解内边界点(i,q,r)，(i,q+1,r)，(i,q,r+1)和(i,q+1,r+1)，i=p+2,p+3,…,l；利用格式A14求解内边界点(p,q,k)，(p,q+1,k)，(p+1,q,k)和(p+1,q+1,k)，k=r+2,r+3,…,n；利用格式A15求解内边界点(p,q,k)，(p,q+1,k)，(p+1,q,k)和(p+1,q+1,k)，k=r-1,r-2,…,1。

利用格式A8～A1依次对子域Ω1～Ω8中的网格点由里向外进行迭代。

3)验证收敛性，若收敛，转到步骤4)，否则回到步骤1)步继续迭代计算。

4)算法停止。

3 数值算例

3.1柱形IEC装置计算

(1)问题描述

为了对该算法的准确性与快速性进行验证，现考虑如下三维泊松问题，以圆柱形惯性静电约束电推力器(IECT)模型为实例进行模拟，不考虑等离子体的影响。

图3 柱形IEC装置电势分布Fig.3 Potential distribution of cylindrical IEC device

(2)结果比较

对该方法的快速性进行研究，将其与现在常用的泊松方程迭代方法进行对比分析，结果如表1所示。

表1 多种迭代格式在不同终止条件下的求解时间Table 1 Solution time of various iterative schemes under different termination conditions ms

注：超松弛因子取为ω=1.6，为模拟求得的较快系数，可以尝试其他值。

从表1可以看出，对于不同的迭代终止条件，P-SOR相较于其他算法的求解时间均大幅领先，可见其快速性是较为稳定的。

3.2 球形IEC装置计算

(1)问题描述

对该算法适用的广泛性开展进一步研究，参考斯图加特大学对于球形IEC装置的研究进展[16]，该装置如图4所示，对该装置进行简化模拟，不考虑等离子体的影响。

图4 球形IEC装置Fig.4 Spherical IEC device

(2)结果分析

使用该P-SOR算法对于柱状或者球状等结构都能较为精确地计算出与实际情况一致的结果，并与串行运算结果相符。

图5 球形IEC装置电势分布Fig.5 Potential distribution of spherical IEC device

进一步地，对超松弛系数ω进行研究，统计多次计算中ω对收敛所需步数的影响，结果如图6所示。可以看出，超松弛系数的选取对求解速度的影响十分显著，和标准形式的超松弛收敛情况类似，极值出现在0<ω<2范围内，本例在取ω值为1.6左右收敛速度最快。而在极值点之后可能会出现不收敛的情况，所以应慎重选取合适的超松弛系数。

图6 迭代步数随超松弛因子的变化关系Fig.6 Relationship between number of iterations and relaxation factor

4 结束语

本文给出的新型P-SOR迭代格式可以很好地对电推力器中复杂电场三维区域进行迭代求解，具有经典SOR格式所没有的显式与并行性。较传统迭代算法快速性显著，至少快出1/2以上，对于求解航天领域研究中的大型三维泊松问题具有非常重要的现实意义。

当然，该方法现阶段尚存在不足，其求解时间距离理论值有一定差距，没有达到只用一种格式求解时间的1/8。究其原因是该并行算法较经典算法复杂许多，其串行计算部分也占据不小的比重，这都导致对其计算时间造成影响。下一步工作需要对算法进行优化，减小串行求解部分占用的比例。而且随着所调用的线程数不断增多，计算机性能越来越好，可以对该方法进行延伸，分解出更多的区域并行求解，使求解速度进一步提高。

References)

[1] TSKHAKAYA D.The particle-in-cell method[J].Contri-butions to Plasma Physics，2010，47(8-9)：563-594.

[2] 陈茂林，夏广庆，魏延明，等.电动帆平行双导线鞘层特性与受力分析[J].物理学报，2016，65(20)：279-289.

CHEN M L，XIA G Q，WEI Y M，et al.Charateristics and stress analysis of sheath of parallel conducting tethers for the electric sail[J].Acta Physica Sinica，2016，65(20)：279-289(in Chinese).

[3] 刘辉，罗晓明，温正，等.GEO卫星霍尔推力器羽流防护结构混合PIC模拟[J].中国空间科学技术，2016，36(1)：63-69.

LIU H，LUO X M，WEN Z，et al.Hybrid-PIC simulation of Hall thruster plume shield on GEO satellites[J].Chinese Space Science and Technology，2016，36(1)：63-69(in Chinese).

[4] 孙安邦，毛根旺，陈茂林，等.霍尔效应推力器羽流的PIC/MCC模拟[J].固体火箭技术，2009，32(6)：638-643.

SUN A B，MAO G W，CHEN M L，et al.PIC/MCC simulation of Hall effect thruster plume[J].Journal of Solid Rocket Technology，2009，32(6)：638-643.

[5] MARTINS S F，FONSECA R A，LU W，et al.Exploring laser-wakefield-accelerator regimes for near-term lasers using particle-in-cell simulation in Lorentz-boosted frames[J].Nature Physics，2010，6(4)： 311.

[6] SAMEH A H.On Jacobi and Jacobi-like algorithms for a parallel computer[J].Mathematics of Computation，1971，25(115)： 579-590.

[7] 蔡放，方逵，李彬.具有完全并行度的SOR算法[J].高等学校计算数学学报，2002，24(1)：11-14.

CAI F，FANG K，LI B.Parallel computation with the perfect degree of parallelism for SOR [J].Numerical Mathematics A Journal of Chinese Universities，2002(1)：11-14.

[8] 岳小宁，肖炳甲，罗正平.基于CUDA的二维泊松方程快速直接求解[J].计算机科学，2013，40(10)：21-23.

YUE X N，XIAO B J，LUO Z P.Fast 2-dimension Poisson direct solver based on CUDA[J].Computer Science，2013，40(10)：21-23(in Chinese).

[9] YOUNG D M，KINCAID D R.A new class of parallel alternating-type iterative methods[J].Journal of Computational & Applied Mathematics，1970，74(1-2)：331-344.

[10] 马月珍，葛永斌，王燕.三维泊松方程基于Richardson外推法的高阶紧致差分方法[J].数学的实践与认识，2011，41(8)：146-152.

MA Y Z，GE Y B，WANG Y.A high-order compact difference method based on Richardson extrapolation for the 3D Poisson equation[J].Mathematics in Practice & Theory，2011，8(8)：146-152(in Chinese).

[11] QUINN M J.Parallel programming in C with MPI and OpenMP[M].New York：McGraw-Hill，2003： 436-513.

[12] MOTOYASU T，NAMBA S，TAKIYAMA K.Mea-surements of localized potential profiles by LIF polarization spectroscopy in an inertial-electrostatic confinement discharge[J].Journal of the Korean Physical Society，2014，65(8)：1205-1208.

[13] 豆桂芳，吴振远，杜艳林.三维泊松方程的七点差分格式[J].工程地球物理学报，2009，6(6)：802-805.

DOU G F，WU Z Y，DU Y L.The seven-point difference scheme for the three-dimensional Poisson's equation[J].Chinese Journal of Engineering Geophysics，2009，6(6)：802-805(in Chinese).

[14] 许秋燕.二维泊松方程和扩散方程的一类显式并行算法[D].济南:山东大学，2010：2-14.

XU Q Y.A class of explicit parallel algorithms for 2d Poisson and diffusion equations[D].Jinan： Shandong University，2010：2-14(in Chinese).

[15] 李春光，徐成贤.确定SOR最佳松弛因子的一个实用算法[J].计算力学学报，2002，19(3)：299-302.

LI C G，XU C X.A practical algorithm determining the optimal relaxation factor of SOR iterative methods[J].Chinese Journal of Computational Mechanics，2002，19(3)：299-302(in Chinese).

[16] CHAN Y A，HERDRICH G，SYRING C，et al.An approach for thrust and losses in inertial electrostatic confinement devices for electric propulsion applications[C].International Electric Propulsion Conference.Kobe：IEPC，2015：083.

(编辑：车晓玲)

Successiveoverrelaxationparallelalgorithmforsolvingcomplexelectricfieldinelectricthrusters

HAN Yajie1，XIA Guangqing1，*，WU Qiuyun1，CHEN Liuwei1，ZHOU Niandong1，WEN Zheng2

1.StateKeyLaboratoryofStructureAnalysisforIndustrialEquipment，DalianUniversityofTechnology，Dalian116024，China2.InstituteofTelecommunicationSatellite，ChinaAcademyofSpaceTechnology，Beijing100094，China

At present，solving the three-dimensional Poisson equation of complex electric field in the electric thruster still has the disadvantages of slowness and poor parallelism.A successive over relaxation parallel (P-SOR)iterative algorithm was proposed to solve the three-dimensional Poisson equation based on the idea of domain decomposition.The solution region was divided into several sub-domains.And a number of explicit schemes were constructed by using the SOR method.Then according to the different schemes，the solving equations for different iterative steps were given for parallel computation.The numerical simulation results show that the calculation time is reduced by at least half.And the P-SOR method has made great progress in the computation speed compared with the traditional iterative scheme.

electric thruster；Poisson equation；successive over relaxation iteration；parallel algorithm

http://zgkj.cast.cn

10.16708/j.cnki.1000-758X.2017.0074

V430

A

2017-03-21；

2017-08-09；录用日期2017-09-12；< class="emphasis_bold">网络出版时间

时间：2017-09-24 16:01:22

http://kns.cnki.net/kcms/detail/11.1859.V.20170924.1601.012.html

国家自然科学基金(11675040)；辽宁省自然科学基金(201602175)；军委科技委H863项目(17-H863-03-ZT-005-072-01)

韩亚杰(1993-)，男，硕士研究生，442300107@qq.com，研究方向为电推进数值模拟与计算

*通讯作者：夏广庆(1979-)，男，副教授，gq.xia@dlut.edu.cn，研究方向为电推进技术

韩亚杰，夏广庆，吴秋云，等.电推力器中复杂电场求解的超松弛并行算法[J].中国空间科学技术，2017，37(5)：68-74.HANYJ，XIAGQ，WUQY，etal.Successiveoverrelaxationparallelalgorithmforsolvingcomplexelectricfieldinelectricthrusters[J].ChineseSpaceScienceandTechnology，2017，37(5)：68-74(inChinese).