Particle-in-Cell等離子體模擬方法的優化研究

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

(蘭州空間技術物理研究所，真空低溫技術與物理重點實驗室，甘肅蘭州730000)

1 引言

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

近年來，PIC方法作為一種重要的模擬方法已經開始應用到航天器與空間等離子體的相互作用中［3］。但是在空間等離子體與航天器相互作用的充放電效應仿真分析中，由于模擬空間大、航天器結構復雜等因素，導致整個計算量非常大，限制了實際工程應用。因此，針對航天器充放電效應的PIC方法仿真研究，研究優化其計算方法，對實際工程應用有重要意義。

2 PIC方法的計算量

2.1 PIC 方法

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

圖1 PIC模擬流程圖

2.2 PIC方法各過程的計算量分析

在PIC方法模擬的每個循環內，一般都包括電荷權重、解Maxwell方程、粒子的運動等三個步驟。本文針對PIC方法的計算量研究，選擇的是二維空間其模擬參數如表1所示。

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

表1 GEO軌道上的磁暴等離子體環境參數

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

c)網格大小為 1×1 mm，時間步長為 Δt=2×10－8s，計算精度設為1×10－5;

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

圖2 模擬空間的邊界設置

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

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

對于Poisson方程，在二維空間內，可展開為:

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

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

為了更加直觀地分析PIC各部分所耗時間的百分比，分別以計算區域為200×200 m和1000×1000 m兩種情況來分析，圖3是根據模擬計算時，記錄的各部分的模擬耗時繪制的柱狀圖。

圖3 PIC模擬各部分所耗時間的柱狀圖

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

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

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

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

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

圖4 計算區域為250×250m時，迭代次數與松弛因子ω的關系

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

對于特定條件，如正方形區域的Dirichlet問題，ω0可使用經驗公式4得到:

式中 l+1為正方形區域中任意一邊的節點數。

如果計算區域為矩形區域，并采用正方形網格分割時，可使用經驗公式5選取。圖4(b)就是當計算區域為正方形時，按公式5畫出的網格數與最佳松弛因子的關系圖。當l=m=250時，利用公式5可以得到ω0=1.975，這與圖4(a)得到的結果一致。

式中 l+1和m+1分別為矩形區域兩邊的節點數。

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

3.2 初始化矩陣對收斂速度的影響

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

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

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

圖5 采用零矩陣作為每個時間步長內的初始矩陣

圖6 采用上一個時間步長內得到的解矩陣作為初始矩陣

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

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

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

同時，當選擇的并非最佳松弛因子時，采用上一步的解矩陣所節約的計算量會更大。這一結論可以通過同樣的模擬計算對比得到。

4 結論

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

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

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

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

［2］王虹宇 王友年等.Particle－in－Cell模擬的發展:物理考慮和計算技術［R］.大連理工大學博士后工作報告，2010:1－3.

［3］買勝利，王立等.材料二次電子發射特性對表面充電過程影響的數值模擬研究［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］王元淳.邊界元法基礎［M］.上海交通大學出版社，1988:1－3.

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