電推力器中復(fù)雜電場(chǎng)求解的超松弛并行算法

，，*，，，，

1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024 2.中國(guó)空間技術(shù)研究院 通信衛(wèi)星事業(yè)部，北京 100094

電推力器中復(fù)雜電場(chǎng)求解的超松弛并行算法

韓亞杰1，夏廣慶1，*，吳秋云1，陳留偉1，周念東1，溫正2

1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024 2.中國(guó)空間技術(shù)研究院 通信衛(wèi)星事業(yè)部，北京 100094

目前，電推力器中復(fù)雜電場(chǎng)的三維泊松方程求解仍存在速度偏慢、并行性欠佳等缺陷。通過(guò)分解三維求解區(qū)域，針對(duì)復(fù)雜電場(chǎng)提出一種超松弛并行(P-SOR)迭代算法。將求解區(qū)域劃分為多個(gè)子塊，利用超松弛迭代格式構(gòu)造出若干分組顯式格式，分別給出不同迭代步數(shù)下的求解方程以進(jìn)行并行計(jì)算。通過(guò)數(shù)值模擬，計(jì)算時(shí)間至少縮短一半以上，該P(yáng)-SOR方法較傳統(tǒng)迭代格式在快速性方面取得較大進(jìn)展，對(duì)電推進(jìn)領(lǐng)域的數(shù)值仿真研究起到促進(jìn)作用。

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

泊松方程是等離子體物理、靜電學(xué)中常見(jiàn)的偏微分方程，其應(yīng)用范圍較廣，特別是在電推力器中復(fù)雜電場(chǎng)的求解方面?，F(xiàn)階段，大規(guī)模的科學(xué)工程計(jì)算對(duì)于速度快、容量大的并行機(jī)和有效的數(shù)值并行算法均有較大的需求，特別是在航天領(lǐng)域。基于粒子云網(wǎng)格(Particle In Cell，PIC)方法[1]的數(shù)值計(jì)算，通過(guò)泊松方程來(lái)求解復(fù)雜電場(chǎng)，在電帆仿真計(jì)算[2]，霍爾推力器的相關(guān)研究[3]等領(lǐng)域均有廣泛應(yīng)用。但由于三維模型求解速度很慢，諸多研究人員針對(duì)模型做出簡(jiǎn)化，如采用二維模型、二維軸對(duì)稱模型[4-5]等，其精確度難免會(huì)下降。所以本文針對(duì)加快三維泊松方程的求解速度開(kāi)展研究。

國(guó)內(nèi)外已有一些專家學(xué)者在該領(lǐng)域進(jìn)行了探索?，F(xiàn)在通用的并行模塊大多采用Jacobi迭代[6]，其具有顯性的并行性。但是對(duì)于求解速度，超松弛迭代等方法則要遠(yuǎn)遠(yuǎn)高于Jacobi格式，也有很多學(xué)者對(duì)其進(jìn)行研究。例如，文獻(xiàn)[7]提出一種針對(duì)點(diǎn)的SOR并行方法，但是在擴(kuò)展性方面存在較大的劣勢(shì)，不利于大規(guī)模應(yīng)用。文獻(xiàn)[8]提出一種新型的、應(yīng)用于統(tǒng)一計(jì)算架構(gòu)(CUDA)的并行算法，利用GPU的硬件特性對(duì)泊松方程的求解進(jìn)行加速。文獻(xiàn)[9]提出一種塊SOR理論，改進(jìn)其求解速度，并對(duì)收斂性做出研究。文獻(xiàn)[10]提出一種基于Richardson外推法的三維泊松方程差分格式，具有使用網(wǎng)格基架點(diǎn)少、精度高、穩(wěn)定性好等優(yōu)點(diǎn)。但是對(duì)于三維泊松方程的多線程并行計(jì)算的研究仍較為欠缺。

本文提出一種基于三維泊松方程的并行計(jì)算理論，適用于OpenMP并行模塊[11],并采用電場(chǎng)結(jié)構(gòu)較為復(fù)雜的柱形和球形慣性靜電約束電推力器[12](Inertial Electrostatic Confinement Thruster，IECT)為模型進(jìn)行相應(yīng)的數(shù)值模擬，通過(guò)對(duì)比不同迭代方法的求解時(shí)間，對(duì)其快速性進(jìn)行驗(yàn)證。

1 泊松問(wèn)題及分組顯式格式

1.1泊松問(wèn)題

考慮在長(zhǎng)方體區(qū)域Ω上的三維泊松方程，

式中：Ω=[0≤x≤L]×[0≤y≤M]×[0≤z≤N]為求解區(qū)域，?Ω為Ω的邊界；u(x,y,z)為待求函數(shù)，f(x,y,z)和g(x,y,z)為已知光滑函數(shù)。

對(duì)于泊松方程內(nèi)部點(diǎn)的逼近方法有很多，這里采用最為常見(jiàn)的七點(diǎn)差分格式[13]：

對(duì)區(qū)域Ω進(jìn)行一致剖分，空間步長(zhǎng)取為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方法很符合求解此類問(wèn)題，但它沒(méi)有明顯的并行性，所以這里將給出一個(gè)新的并行SOR迭代算法，簡(jiǎn)稱為P-SOR。根據(jù)對(duì)二維并行迭代格式的研究[14]，推導(dǎo)出基于三維泊松方程的8個(gè)基本SOR迭代格式：

式中：ω為松弛因子。

1.2 分組顯式格式

首先，考慮求解域中相鄰8個(gè)網(wǎng)格點(diǎn)的數(shù)值解，如圖1所示，聯(lián)立A1～A8可形成(8×8)階方程組，采用分組顯式的思想，即可得到求解域節(jié)點(diǎn)(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)處的數(shù)值解公式。

圖1 中間點(diǎn)迭代示意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ω

類似地，利用分組顯式基于“四點(diǎn)一組”的思想可構(gòu)造其他6個(gè)新的格式。

利用格式A1,A3,A5,A7求解內(nèi)邊界點(diǎn)(i,j0,k)，(i+1,j0,k)，(i,j0,k+1)和(i+1,j0,k+1)，j0=j+2,j+3,…,m可得格式A10：

用類似的方法即可求出另外5個(gè)迭代格式A11～A15，其基本表現(xiàn)形式相同，這里不再詳述。

特別地，當(dāng)f(i,j,k)=0時(shí)，該泊松方程即退化為拉普拉斯方程。對(duì)于超松弛系數(shù)ω的選取，可以參考文獻(xiàn)[15]。

2 P-SOR算法流程

2.1區(qū)域分解

首先，將區(qū)域Ω分解為若干子域，這里以8個(gè)非重疊子域?yàn)槔?其他以此類推)，如圖2所示，用6個(gè)網(wǎng)格平面(i=p,p+1,j=q,q+1,k=r,r+1)將求解區(qū)域劃分的8個(gè)子域，其中p,q,r為任意正整數(shù)，0

圖2 求解域網(wǎng)格劃分Fig.2 Meshing partitions in the calculation domain

2.2 迭代算法

顯然，在求解過(guò)程中對(duì)于相同的區(qū)域不能使用同一種迭代格式，這樣會(huì)造成不同時(shí)間步的電勢(shì)值無(wú)法更新迭代。所以，對(duì)于奇數(shù)次迭代和偶數(shù)次迭代，采用不同的迭代算法進(jìn)行實(shí)現(xiàn)，步驟如下：

1)迭代次數(shù)為奇數(shù)。

利用格式A1～A8依次對(duì)子域Ω1～Ω8中的網(wǎng)格點(diǎn)由外向內(nèi)進(jìn)行迭代。

2)迭代次數(shù)為偶數(shù)。

利用格式A9求解網(wǎng)格中心交叉點(diǎn)(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求解內(nèi)邊界點(diǎn)(p,j,r)，(p+1,j,r)，(p,j,r+1)和(p+1,j,r+1)，j=q+2,q+3,…,m；利用格式A11求解內(nèi)邊界點(diǎn)(p,j,r)，(p+1,j,r)，(p,j,r+1)和(p+1,j,r+1)，j=q-1,q-2,…,1；利用格式A12求解內(nèi)邊界點(diǎn)(i,q,r)，(i,q+1,r)，(i,q,r+1)和(i,q+1,r+1)，i=p-1,p-2,…,1；利用格式A13求解內(nèi)邊界點(diǎn)(i,q,r)，(i,q+1,r)，(i,q,r+1)和(i,q+1,r+1)，i=p+2,p+3,…,l；利用格式A14求解內(nèi)邊界點(diǎn)(p,q,k)，(p,q+1,k)，(p+1,q,k)和(p+1,q+1,k)，k=r+2,r+3,…,n；利用格式A15求解內(nèi)邊界點(diǎn)(p,q,k)，(p,q+1,k)，(p+1,q,k)和(p+1,q+1,k)，k=r-1,r-2,…,1。

利用格式A8～A1依次對(duì)子域Ω1～Ω8中的網(wǎng)格點(diǎn)由里向外進(jìn)行迭代。

3)驗(yàn)證收斂性，若收斂，轉(zhuǎn)到步驟4)，否則回到步驟1)步繼續(xù)迭代計(jì)算。

4)算法停止。

3 數(shù)值算例

3.1柱形IEC裝置計(jì)算

(1)問(wèn)題描述

為了對(duì)該算法的準(zhǔn)確性與快速性進(jìn)行驗(yàn)證，現(xiàn)考慮如下三維泊松問(wèn)題，以圓柱形慣性靜電約束電推力器(IECT)模型為實(shí)例進(jìn)行模擬，不考慮等離子體的影響。

圖3 柱形IEC裝置電勢(shì)分布Fig.3 Potential distribution of cylindrical IEC device

(2)結(jié)果比較

對(duì)該方法的快速性進(jìn)行研究，將其與現(xiàn)在常用的泊松方程迭代方法進(jìn)行對(duì)比分析，結(jié)果如表1所示。

表1 多種迭代格式在不同終止條件下的求解時(shí)間Table 1 Solution time of various iterative schemes under different termination conditions ms

注：超松弛因子取為ω=1.6，為模擬求得的較快系數(shù)，可以嘗試其他值。

從表1可以看出，對(duì)于不同的迭代終止條件，P-SOR相較于其他算法的求解時(shí)間均大幅領(lǐng)先，可見(jiàn)其快速性是較為穩(wěn)定的。

3.2 球形IEC裝置計(jì)算

(1)問(wèn)題描述

對(duì)該算法適用的廣泛性開(kāi)展進(jìn)一步研究，參考斯圖加特大學(xué)對(duì)于球形IEC裝置的研究進(jìn)展[16]，該裝置如圖4所示，對(duì)該裝置進(jìn)行簡(jiǎn)化模擬，不考慮等離子體的影響。

圖4 球形IEC裝置Fig.4 Spherical IEC device

(2)結(jié)果分析

使用該P(yáng)-SOR算法對(duì)于柱狀或者球狀等結(jié)構(gòu)都能較為精確地計(jì)算出與實(shí)際情況一致的結(jié)果，并與串行運(yùn)算結(jié)果相符。

圖5 球形IEC裝置電勢(shì)分布Fig.5 Potential distribution of spherical IEC device

進(jìn)一步地，對(duì)超松弛系數(shù)ω進(jìn)行研究，統(tǒng)計(jì)多次計(jì)算中ω對(duì)收斂所需步數(shù)的影響，結(jié)果如圖6所示?？梢钥闯觯沙谙禂?shù)的選取對(duì)求解速度的影響十分顯著，和標(biāo)準(zhǔn)形式的超松弛收斂情況類似，極值出現(xiàn)在0<ω<2范圍內(nèi)，本例在取ω值為1.6左右收斂速度最快。而在極值點(diǎn)之后可能會(huì)出現(xiàn)不收斂的情況，所以應(yīng)慎重選取合適的超松弛系數(shù)。

圖6 迭代步數(shù)隨超松弛因子的變化關(guān)系Fig.6 Relationship between number of iterations and relaxation factor

4 結(jié)束語(yǔ)

本文給出的新型P-SOR迭代格式可以很好地對(duì)電推力器中復(fù)雜電場(chǎng)三維區(qū)域進(jìn)行迭代求解，具有經(jīng)典SOR格式所沒(méi)有的顯式與并行性。較傳統(tǒng)迭代算法快速性顯著，至少快出1/2以上，對(duì)于求解航天領(lǐng)域研究中的大型三維泊松問(wèn)題具有非常重要的現(xiàn)實(shí)意義。

當(dāng)然，該方法現(xiàn)階段尚存在不足，其求解時(shí)間距離理論值有一定差距，沒(méi)有達(dá)到只用一種格式求解時(shí)間的1/8。究其原因是該并行算法較經(jīng)典算法復(fù)雜許多，其串行計(jì)算部分也占據(jù)不小的比重，這都導(dǎo)致對(duì)其計(jì)算時(shí)間造成影響。下一步工作需要對(duì)算法進(jìn)行優(yōu)化，減小串行求解部分占用的比例。而且隨著所調(diào)用的線程數(shù)不斷增多，計(jì)算機(jī)性能越來(lái)越好，可以對(duì)該方法進(jìn)行延伸，分解出更多的區(qū)域并行求解，使求解速度進(jìn)一步提高。

References)

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

[2] 陳茂林，夏廣慶，魏延明，等.電動(dòng)帆平行雙導(dǎo)線鞘層特性與受力分析[J].物理學(xué)報(bào)，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衛(wèi)星霍爾推力器羽流防護(hù)結(jié)構(gòu)混合PIC模擬[J].中國(guó)空間科學(xué)技術(shù)，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] 孫安邦，毛根旺，陳茂林，等.霍爾效應(yīng)推力器羽流的PIC/MCC模擬[J].固體火箭技術(shù)，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，F(xiàn)ONSECA 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].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，2002，24(1)：11-14.

CAI F，F(xiàn)ANG 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].計(jì)算機(jī)科學(xué)，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].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，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] 豆桂芳，吳振遠(yuǎn)，杜艷林.三維泊松方程的七點(diǎn)差分格式[J].工程地球物理學(xué)報(bào)，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] 許秋燕.二維泊松方程和擴(kuò)散方程的一類顯式并行算法[D].濟(jì)南:山東大學(xué)，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最佳松弛因子的一個(gè)實(shí)用算法[J].計(jì)算力學(xué)學(xué)報(bào)，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">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2017-09-24 16:01:22

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

國(guó)家自然科學(xué)基金(11675040)；遼寧省自然科學(xué)基金(201602175)；軍委科技委H863項(xiàng)目(17-H863-03-ZT-005-072-01)

韓亞杰(1993-)，男，碩士研究生，442300107@qq.com，研究方向?yàn)殡娡七M(jìn)數(shù)值模擬與計(jì)算

*通訊作者：夏廣慶(1979-)，男，副教授，gq.xia@dlut.edu.cn，研究方向?yàn)殡娡七M(jìn)技術(shù)

韓亞杰，夏廣慶，吳秋云，等.電推力器中復(fù)雜電場(chǎng)求解的超松弛并行算法[J].中國(guó)空間科學(xué)技術(shù)，2017，37(5)：68-74.HANYJ，XIAGQ，WUQY，etal.Successiveoverrelaxationparallelalgorithmforsolvingcomplexelectricfieldinelectricthrusters[J].ChineseSpaceScienceandTechnology，2017，37(5)：68-74(inChinese).