Particle-in-Cell等離子體模擬方法的優(yōu)化研究

劉明正，楊生勝，湯道坦，秦曉剛，顏則東

(蘭州空間技術(shù)物理研究所，真空低溫技術(shù)與物理重點(diǎn)實(shí)驗(yàn)室，甘肅蘭州730000)

1 引言

PIC(Particle-in-Cell)方法［1］是目前研究等離子體集體動(dòng)力學(xué)行為的一種強(qiáng)有力的數(shù)值手段，其通過跟蹤大量帶電粒子在外加及自洽電磁場(chǎng)中的運(yùn)動(dòng)，通過統(tǒng)計(jì)平均得到宏觀特性及運(yùn)動(dòng)規(guī)律的一種方法。從理論上說，只要粒子數(shù)和網(wǎng)格數(shù)量足夠多，PIC方法就可以處理幾乎所有的等離子體問題［2］。該方法廣泛應(yīng)用于等離子體物理所涉及的許多領(lǐng)域，如空間等離子體、自由電子激光、激光等離子體相互作用等。

近年來，PIC方法作為一種重要的模擬方法已經(jīng)開始應(yīng)用到航天器與空間等離子體的相互作用中［3］。但是在空間等離子體與航天器相互作用的充放電效應(yīng)仿真分析中，由于模擬空間大、航天器結(jié)構(gòu)復(fù)雜等因素，導(dǎo)致整個(gè)計(jì)算量非常大，限制了實(shí)際工程應(yīng)用。因此，針對(duì)航天器充放電效應(yīng)的PIC方法仿真研究，研究?jī)?yōu)化其計(jì)算方法，對(duì)實(shí)際工程應(yīng)用有重要意義。

2 PIC方法的計(jì)算量

2.1 PIC 方法

PIC方法模擬等離子體的基本過程是:根據(jù)粒子的相空間(珒r(jià)，珒v)計(jì)算出空間網(wǎng)格上的電流及電荷密度，求解Maxwell方程計(jì)算出空間網(wǎng)格上的電磁場(chǎng)珗E，珗B，使用權(quán)重分配方法并利用洛倫茲方程求出粒子新的相空間，如此形成一個(gè)不斷迭代的循環(huán)。整個(gè)PIC方法的模擬流程如圖1所示。

圖1 PIC模擬流程圖

2.2 PIC方法各過程的計(jì)算量分析

在PIC方法模擬的每個(gè)循環(huán)內(nèi)，一般都包括電荷權(quán)重、解Maxwell方程、粒子的運(yùn)動(dòng)等三個(gè)步驟。本文針對(duì)PIC方法的計(jì)算量研究，選擇的是二維空間其模擬參數(shù)如表1所示。

a)空間環(huán)境為GEO軌道上的磁暴等離子體環(huán)境，具體的等離子體參數(shù)見表1［4］。

表1 GEO軌道上的磁暴等離子體環(huán)境參數(shù)

b)為了方便研究，模擬空間區(qū)域和衛(wèi)星形狀均采用正方形，且衛(wèi)星尺寸是10×10m，位于整個(gè)模擬區(qū)域的中心;

c)網(wǎng)格大小為 1×1 mm，時(shí)間步長(zhǎng)為 Δt=2×10－8s，計(jì)算精度設(shè)為1×10－5;

d)邊界條件:模擬區(qū)域的外邊界電勢(shì)為零;衛(wèi)星周圍邊界條件是:%n=－ρs/ε0，如圖2所示。

圖2 模擬空間的邊界設(shè)置

由于GEO軌道環(huán)境，磁場(chǎng)強(qiáng)度很小，可以忽略。所以一般采用PIC的靜電模型來模擬GEO軌道上的等離子環(huán)境。這時(shí)，空間中的Maxwell方程可簡(jiǎn)化為Poisson方程:

解Poisson方程一般有數(shù)值差分法(Finite Difference Method，簡(jiǎn)稱FDM)、有限元法(Finite Element Method，簡(jiǎn)稱 FEM)和邊界元法(Boundary Element Method，簡(jiǎn)稱BEM)三種數(shù)值方法［5］?？紤]到這三種方法的優(yōu)劣，一般選用數(shù)值差分法。

對(duì)于Poisson方程，在二維空間內(nèi)，可展開為:

采用逐次超松弛迭代法(Successive Over Relaxation Method，簡(jiǎn)稱SOR法)解式2的迭代公式為:

式3中的ω為松弛因子，Δx=Δy=h。當(dāng)ω=1時(shí)，上式就是高斯－賽德爾迭代公式［6］。

為了更加直觀地分析PIC各部分所耗時(shí)間的百分比，分別以計(jì)算區(qū)域?yàn)?00×200 m和1000×1000 m兩種情況來分析，圖3是根據(jù)模擬計(jì)算時(shí)，記錄的各部分的模擬耗時(shí)繪制的柱狀圖。

圖3 PIC模擬各部分所耗時(shí)間的柱狀圖

從圖3可以看出，在PIC模擬的每個(gè)循環(huán)內(nèi)，計(jì)算量主要集中在求解Poisson方程這一步。當(dāng)計(jì)算區(qū)域尺寸為200×200 m時(shí)，求解Poisson方程大約占用了總時(shí)間的91%;而當(dāng)計(jì)算區(qū)域尺寸為1000×1000 m時(shí)，該部分幾乎占用了全部的計(jì)算時(shí)間。由此可見，一個(gè)好的、高效的Poisson方程求解算法可以顯著地降低PIC模擬的計(jì)算量、縮短計(jì)算時(shí)間。

3 影響Poisson方程求解速度的因素

當(dāng)其他條件不變時(shí)，SOR法解Poisson方程的收斂速度取決于所選擇的松弛因子和初始化的電勢(shì)矩陣兩個(gè)因素。下面著重討論這兩個(gè)因素對(duì)計(jì)算速度的影響。

3.1 松弛因子對(duì)收斂速度的影響及選取方法

為了保證SOR方法收斂，松弛因子ω必須滿足0＜ω＜2［6］，且通常當(dāng)1＜ω＜2時(shí)，收斂速度較快。為了分析ω對(duì)收斂速度的影響，這里選擇計(jì)算區(qū)域?yàn)?50×250 m的模擬區(qū)域尺寸作為研究對(duì)象，其他等離子體參數(shù)不變。圖4(a)就是在上述條件下得到的松弛因子與解Poisson方程所需迭代次數(shù)的關(guān)系圖。

圖4 計(jì)算區(qū)域?yàn)?50×250m時(shí)，迭代次數(shù)與松弛因子ω的關(guān)系

從圖4(a)可以看出，松弛因子對(duì)求解速度的影響很大:當(dāng)ω=1時(shí)，得到Poisson方程的數(shù)值解需要24800個(gè)迭代循環(huán)，而當(dāng)ω=1.975時(shí)，得到該結(jié)果只需502步，兩者相差近50倍。

對(duì)于特定條件，如正方形區(qū)域的Dirichlet問題，ω0可使用經(jīng)驗(yàn)公式4得到:

式中 l+1為正方形區(qū)域中任意一邊的節(jié)點(diǎn)數(shù)。

如果計(jì)算區(qū)域?yàn)榫匦螀^(qū)域，并采用正方形網(wǎng)格分割時(shí)，可使用經(jīng)驗(yàn)公式5選取。圖4(b)就是當(dāng)計(jì)算區(qū)域?yàn)檎叫螘r(shí)，按公式5畫出的網(wǎng)格數(shù)與最佳松弛因子的關(guān)系圖。當(dāng)l=m=250時(shí)，利用公式5可以得到ω0=1.975，這與圖4(a)得到的結(jié)果一致。

式中 l+1和m+1分別為矩形區(qū)域兩邊的節(jié)點(diǎn)數(shù)。

式4和式5提供了在指定條件下得到最佳松弛因子的經(jīng)驗(yàn)公式。然而并不符合這些特定條件，最佳松弛因子ω0只能靠經(jīng)驗(yàn)選取。對(duì)比圖4(a)和圖4(b)可以發(fā)現(xiàn)，采用枚舉方法選取的ω0與公式獲得的結(jié)果相同，因此對(duì)于一般情況，可以通過在第一個(gè)時(shí)間步長(zhǎng)內(nèi)，不同ω0下求解Poisson方程的速度對(duì)比分析，選出整個(gè)PIC模擬過程中的最佳松弛因子。

3.2 初始化矩陣對(duì)收斂速度的影響

使用SOR法解差分方程時(shí)，還需要一個(gè)初始矩陣。從理論上來說，初始矩陣與解矩陣越接近，需要的迭代次數(shù)越少。但是，由于在得到數(shù)值解之前，一般無法估計(jì)解矩陣，所以常選擇零矩陣作為初始矩陣。對(duì)于一些具體問題，需要多次求解Poisson方程，如果始終以零矩陣作為初始矩陣，雖然仍能得到正確結(jié)果，但是計(jì)算量并非最小。對(duì)應(yīng)地，需要的計(jì)算時(shí)間也比較長(zhǎng)。

在PIC方法中，每一個(gè)時(shí)間步長(zhǎng)內(nèi)都需要根據(jù)當(dāng)前空間中粒子的相空間計(jì)算該區(qū)域中的電勢(shì)。而為了保證PIC方法的數(shù)值穩(wěn)定性，要求粒子在一個(gè)時(shí)間步長(zhǎng)內(nèi)移動(dòng)的距離不大于網(wǎng)格尺寸，這就使得空間中粒子在經(jīng)過一個(gè)時(shí)間步長(zhǎng)后相空間變化較小。所以每一個(gè)時(shí)間步長(zhǎng)內(nèi)空間中的電勢(shì)變化也比較小。這就是說，上一個(gè)時(shí)間步長(zhǎng)時(shí)的空間電勢(shì)與當(dāng)前時(shí)間步長(zhǎng)時(shí)的空間電勢(shì)比較接近。如果每次解Poisson方程時(shí)，都是以上一個(gè)時(shí)間步長(zhǎng)時(shí)的電勢(shì)作為初始矩陣，理論上可以大大減少計(jì)算量。

這里同樣采用PIC方法模擬上文提到的等離子體環(huán)境和網(wǎng)格等參數(shù)，計(jì)算區(qū)域尺寸為250×250m。為了便于比較，這里將松弛因子ω定為1.975。如果選擇零矩陣作為初始矩陣，則在每一個(gè)時(shí)間步長(zhǎng)內(nèi)，解Poisson方程所用的迭代次數(shù)如圖5(a)所示。對(duì)于的航天器表面介質(zhì)上的最小電勢(shì)隨時(shí)間的變化如圖5(b)所示。如果采用上個(gè)時(shí)間步長(zhǎng)時(shí)得到的解矩陣作為初始矩陣，則PIC方法的每個(gè)循環(huán)內(nèi)，求解Poisson方程所用的迭代次數(shù)和介質(zhì)表面的最小電勢(shì)隨時(shí)間的變化曲線如圖6所示。

圖5 采用零矩陣作為每個(gè)時(shí)間步長(zhǎng)內(nèi)的初始矩陣

圖6 采用上一個(gè)時(shí)間步長(zhǎng)內(nèi)得到的解矩陣作為初始矩陣

從圖5可以看出:如果采用零矩陣作為初始矩陣，迭代次數(shù)隨著模擬時(shí)間逐漸增加，也就是說，計(jì)算量隨著模擬時(shí)間逐漸增加。同時(shí)，介質(zhì)表面的最低電勢(shì)也隨著時(shí)間逐漸降低。對(duì)比圖5(a)和(b)可以發(fā)現(xiàn)，迭代次數(shù)隨最低電勢(shì)的降低而降低。其原因是初始矩陣與解矩陣的差別越來越多，計(jì)算量也越來越大。從圖6可以看出:如果采用上一個(gè)時(shí)間步長(zhǎng)內(nèi)的解矩陣作為SOR法解Poisson方程的初始矩陣，則每次解Poisson方程的計(jì)算量隨著時(shí)間逐漸減少，并且最終穩(wěn)定在某一個(gè)較小的值附近。從圖6(a)中同樣可以發(fā)現(xiàn):在第一個(gè)時(shí)間步長(zhǎng)內(nèi)，迭代步數(shù)很大，這是由于在該時(shí)間步長(zhǎng)內(nèi)采用零矩陣作為初始矩陣，而該初始矩陣一般與解矩陣有較大差別(這里Poisson方程的解矩陣中的最小值為－69.12)，所以所需的迭代次數(shù)較多。

對(duì)比圖5(a)和6(a)可以看出:采用上一步的解矩陣作為初始矩陣比使用零矩陣作為初始矩陣可以顯著降低計(jì)算量。這里大約能節(jié)約近一半的計(jì)算量。

對(duì)比圖5(b)和6(b)可以看出:采用不同的初始條件，得到的介質(zhì)表面最低電勢(shì)隨模擬時(shí)間的變化曲線基本相同，因此初始矩陣的選擇對(duì)最終的結(jié)果沒有影響。

同時(shí)，當(dāng)選擇的并非最佳松弛因子時(shí)，采用上一步的解矩陣所節(jié)約的計(jì)算量會(huì)更大。這一結(jié)論可以通過同樣的模擬計(jì)算對(duì)比得到。

4 結(jié)論

在利用PIC方法模擬航天器與空間等離子體相互作用時(shí)，求解泊松方程是計(jì)算量最大的一個(gè)步驟。因此，提高求解泊松方程的計(jì)算速度，可大大優(yōu)化PIC方法在航天器與空間等離子體相互作用數(shù)值模擬中的應(yīng)用。通過本文的分析研究，針對(duì)提高求解速度得出以下結(jié)論:

(1)選擇合適的松弛因子，可以顯著減少解泊松方程的計(jì)算量。本文給出了一種適用于一般情況的最佳松弛因子選擇方法;

(2)當(dāng)求解泊松方程時(shí)，采用上一個(gè)時(shí)間步長(zhǎng)內(nèi)的解矩陣作為初始矩陣的方法能顯著降低計(jì)算量、提高計(jì)算效率，并且隨著模擬時(shí)間的增長(zhǎng)，效果越來越明顯。

［1］Birdsall C K，Landon A B.Plasma Physics via Computer Simulation［R］.Longdon:Institute of Physics Publishing，1991.

［2］王虹宇 王友年等.Particle－in－Cell模擬的發(fā)展:物理考慮和計(jì)算技術(shù)［R］.大連理工大學(xué)博士后工作報(bào)告，2010:1－3.

［3］買勝利，王立等.材料二次電子發(fā)射特性對(duì)表面充電過程影響的數(shù)值模擬研究［J］.真空與低溫，2006，12(2):87－90.

［4］Anderson B.Jeffrey.Natural Orbital Environment Guidelines for Use in Aerospace Vehicle Development［M］.NASA Marshall Space Flight Center.1994:31－33.

［5］王元淳.邊界元法基礎(chǔ)［M］.上海交通大學(xué)出版社，1988:1－3.

［6］顏慶津.數(shù)值分析(第三版)［M］.北京航空航天大學(xué)出版社2006.