基于化學物質(zhì)釋放的五點傅里葉變換求解泊松方程的方法

徐朝輝 趙海生 許正文 吳健 馮杰 徐彬 馬征征 黃文龍

(中國電波傳播研究所 電波環(huán)境特性及模化技術(shù)重點實驗室,青島 266107)

引 言

在電離層化學物質(zhì)釋放引起的等離子體小尺度不均勻體數(shù)值模擬中,電勢求解是非常重要的一環(huán)．在靜電粒子模型中,靜電力起著決定作用,已知電荷密度分布的情況下,不必求解復雜的麥克斯韋方程組,通過求解泊松方程,就能夠獲得空間電勢分布,進而得到空間電場分布,可用于描述化學物質(zhì)釋放產(chǎn)生等離子體不均勻結(jié)構(gòu)的動力學過程．通常求解泊松方程的方法有矩陣法、迭代法和譜方法[1-2]．

等離子體主要表現(xiàn)為集體運動特性,在等離子體粒子模擬中使用有限大小粒子可以使近距離碰撞作用極大降低,而不改變遠程相互作用,從而使集體運動特性保留下來．粒子模擬中有限大小粒子的電荷不再是集中一點,而是按照一定的形狀連續(xù)地分布在它的有限空間里,因而力的計算方法有很大的不同．因為力的計算采用卷積的形式,顯然沒有必要再像點粒子那樣用兩粒子作用力疊加的辦法來計算,而用傅里葉變換的方法計算要簡單很多,也就是說,求解泊松方程,不必再去使用精確的電荷密度,而只要將電荷密度在空間網(wǎng)格點上作近似的展開即可[3]．這樣的展開一方面收斂快,另一方面不必使用偶極展開項,并且傅里葉變換方法編程程序簡單,占用計算機處理器內(nèi)存少,節(jié)約了粒子模擬的運算時間[4]．

快速傅里葉變換(Fast Fourier Transform,FFT)在空間網(wǎng)格點上求解泊松方程最早由Cooley、Tukey提出[5],后來Weaver在Cooley的基礎(chǔ)上做了詳細的離散傅里葉分析[6]．Birdsall 和 Langdon提出了經(jīng)典的FFT解泊松方程的方法[7],這個方法推導方便,計算效率高,程序?qū)崿F(xiàn)簡單,成為國內(nèi)外沿用至今的方法．近年來,隨著離子模擬應用場景的推廣,根據(jù)各種應用條件的不同也發(fā)展了幾種其他的FFT求解泊松方程的方法[2]．

在將離子模擬方法引入電離層化學物質(zhì)釋放的數(shù)值模擬過程中,我們發(fā)現(xiàn)Birdsall方法在計算精度上有待提高,參考國內(nèi)外相關(guān)文獻[1,5-6],在經(jīng)典譜方法的基礎(chǔ)上[8],提出了一種新的FFT求解泊松方程的方法．結(jié)果表明,改進的方法精度更高,電勢分布更符合解析值,為等離子體粒子模擬過程中靜電模型的泊松方程求解提供了一種新的方法．

1 求解原理

首先給出電離層化學物質(zhì)釋放后系統(tǒng)中靜電場電勢的泊松方程:

(1)

ρ=ρe+ρx+ρO.

(2)

式中:ρ是電子和離子總的電荷密度;下標e、x、O分別代表電子、負離子和氧離子;ε0是真空介電常數(shù)．為了便于計算,這里模擬的系統(tǒng)是釋放早期,約等于或少于負離子回旋周期的時間內(nèi)電子和離子行為,但是忽略了擴散效應和電子吸附過程等其他因素的影響,并非是在整個釋放穩(wěn)定后的狀態(tài)．

1.1 級數(shù)求導方法

級數(shù)求導法是一種常用的泊松方程求解方法[9],該方法的基本思路是通過FFT求得電荷密度和電勢傅里葉級數(shù)的系數(shù),得到電荷密度和電勢的傅里葉級數(shù);將得到的傅里葉級數(shù)代入泊松方程,可得到電勢的傅里葉象函數(shù)對應關(guān)系;對電勢的象函數(shù)做逆傅里葉變換,即可求得電勢．以二維形式為例,對于周期系統(tǒng),在滿足Dirichlet條件下得到電勢的傅里葉象函數(shù):

(3)

式中:Pkx,ky為電荷密度的象函數(shù);Nx和Ny為網(wǎng)格大小．可以看到,電勢象函數(shù)在頻域下為顯式格式,求解十分方便．不過還需要注意到該方法可能引起吉布斯現(xiàn)象[10],所以需要一個額外的濾波器通過減小高階傅里葉系數(shù)的方法來平滑邊界,同時還要注意到對模型化簡時,用到了近似替代的簡化手段,這些操作使計算精度有所下降．

1.2 Birdsall 和 Langdon經(jīng)典解法

Birdsall和Langdon給出了一種經(jīng)典FFT求解泊松方程的方法[5],以二維形式為例,對于周期系統(tǒng),求解泊松方程的基本流程[3]如圖1所示．

圖1 FFT方法求解電勢的流程圖

根據(jù)文獻[5]可以得到,電勢象函數(shù)的表達式為

(4)

可以看到,Birdsall方法中,電勢象函數(shù)的計算過程有了很大簡化,但同時要考慮網(wǎng)格大小Δx和Δy對近似替代精度的影響．

1.3 五點傅里葉變換方法

傳統(tǒng)等離子體粒子模擬中FFT求解泊松方程的方法已經(jīng)十分成熟[11],計算精確度在大多數(shù)情況下能夠滿足粒子模擬的需要．然而在開展電離層化學物質(zhì)釋放產(chǎn)生的電子密度空洞和等離子體云的小尺度擾動數(shù)值模擬仿真過程中,為了進一步精確模擬電子密度空洞邊界不規(guī)則體的演化,得到不穩(wěn)定性相干結(jié)構(gòu),本文提出一種改進的FFT求解泊松方程的方法,滿足了化學物質(zhì)釋放早期演化過程模擬精度的需要．

對于周期邊界條件下的二維系統(tǒng),將式(1)差分得到

(5)

將式(5)用中心差分法繼續(xù)分解,有

(6)

整理得到

(7)

(8)

對式(8)化簡,合并為三角函數(shù)形式有

(9)

在粒子模擬網(wǎng)格劃分時,網(wǎng)格長度設為Δx=Δy,式(9)進一步化簡為

(10)

此時,考慮到對二階偏導做中心差分存在較大誤差,為了彌補二階差分對差分精度的影響我們引入差分誤差修正項Rn(Rn≤ε),

(11)

ε取值根據(jù)網(wǎng)格長度的不同而變化[4, 12]．則電勢的傅里葉象函數(shù)可表示為

(12)

最后,對電勢的傅里葉象函數(shù)求傅里葉逆變換得到電勢φx,y的值．

改進的五點傅里葉變換方法,計算效率方面,只需一次插值、一次FFT和IFFT計算,僅比Birdsall方法和級數(shù)求導法多了一次插值運算．但是在傅里葉空間的系數(shù)運算中并未使用近似替代等處理方法,減小了非線性誤差,同時對二階差分誤差采取了誤差彌補措施,進一步提高了差分精度,同時保證了算法的穩(wěn)定性．這種算法精度高、穩(wěn)定性好,而且計算效率尚可,具有重要的理論意義和實際應用價值．

2 求解對比分析

2.1 單點電荷電勢分布舉例

下面使用三種方法求解一個簡單例子,對比分析五點傅里葉改進方法和另外兩種方法的優(yōu)缺點．首先,建立一個64×64二維網(wǎng)格,在化學物質(zhì)釋放數(shù)值模擬領(lǐng)域,網(wǎng)格邊界取為電子德拜長度,邊界條件采用周期邊界．電子和離子的總電荷密度分布已知,為了便于觀察對比求解結(jié)果,這里用一個網(wǎng)格點(33,33)上20個單位正電荷的數(shù)密度給出,密度值為6.6942×10-13C/m2,如圖2所示．

圖2 點電荷密度二維分布

根據(jù)解析公式,可得到點電荷的電勢

(13)

式中:q為點電荷的帶電量,取值為20個單位正電荷;r為二維空間任一點距點電荷的距離．q>0時,電場中各點的電勢都是正值,且隨r的增加而減小．由式(13)可以得到在網(wǎng)格大小為64×64的系統(tǒng)中電勢的解析分布,如圖3所示,電勢的最大量級在10-5．

圖3 解析方法得到的點電荷電勢二維分布

根據(jù)點電荷密度分布,分別使用三種方法求解泊松方程,得到三種方法下電勢的分布,如圖4、5、6所示．五點傅里葉方法求得的電勢最大值在10-7量級,Birds all方法和級數(shù)求導法求得的電勢最大值在10-2量級,單純從最大值誤差精度看,五點傅里葉方法相比其他兩種方法誤差更小,更加趨近真實值．由圖4可見,點電荷的電勢梯度分布均勻,電勢等值線輪廓近似圓形,電勢關(guān)于點電荷對稱分布,近似符合數(shù)值解析結(jié)果．與圖3結(jié)果相比較,二者的電勢空間分布結(jié)構(gòu)基本一致,電勢值與解析值之間的均方誤差(Mean Square Error, MSE)為1.0539×10-12,電勢的空間分布符合實際,其優(yōu)勢足以彌補小量級精度誤差的缺陷．

圖4 五點傅里葉方法得到的點電荷二維電勢分布

圖5為Birdsall經(jīng)典方法得到的近似電勢分布．電勢分布近似為y=x對稱分布,電勢最大值約為10-2量級,與理論解析分布存在較大誤差,電勢值與解析值之間的MSE達到 1.2378×10-6,同時由于折疊計算導致電勢空間分布失真,與解析電勢分布差異較大．

圖5 Birdsall經(jīng)典解法得到的點電荷電勢分布

由圖6可以看到,級數(shù)求導方法得到的電勢分布存在明顯數(shù)值誤差,電勢最大值趨近10-2量級,而MSE達到8.4134×10-6,電勢空間分布出現(xiàn)負值,與理論解析分布亦存在較大差異．

圖6 級數(shù)求導方法得到的點電荷電勢分布

為了研究三種方法數(shù)值解的精度,我們做了三種方法的誤差對數(shù)(lg(e))曲線,如圖7所示．從圖7可以看到,五點傅里葉方法的誤差對數(shù)遠小于其他兩種方法,且數(shù)值精度有了很大提高．

圖7 三種方法的誤差對數(shù)曲線

2.2 求解化學物質(zhì)釋放實際應用

將五點傅里葉變換方法代入中性氣體化學物質(zhì)釋放產(chǎn)生密度空洞的靜電全粒子模擬例子中,設網(wǎng)格為128×128,網(wǎng)格大小為電子德拜半徑λDe,時間尺度為1/ωpe;周期邊界條件下,電子、離子數(shù)密度分布已知,電子密度為中間區(qū)域為零,四周均勻分布,負離子密度為中間盤狀區(qū)域均勻分布,正離子密度在全區(qū)域均勻分布[8]．在模擬系統(tǒng)啟動時,系統(tǒng)呈整體電中性,電勢能為零,如圖8所示．

圖9為利用五點傅里葉變換方法求得的化學物質(zhì)釋放后電勢分布．在系統(tǒng)運行到ωpet=100時,空洞邊界層開始出現(xiàn)明顯的電勢分界,此時負離子運動到邊界層外,而電子由于磁場的束縛不能進入邊界層,導致負電荷在邊界層外積聚形成內(nèi)外電勢差,從而形成一個外向的電場．從圖9可以看到,電勢分布基本聚集在電子空洞的邊界層,里面電勢高外面電勢低,電勢分布細節(jié)清晰,符合理論推導．進而證明,我們采用五點傅里葉方法求解泊松方程滿足系統(tǒng)的要求．

圖8 初始時刻電子、負離子和氧離子密度分布(順序從上到下)

圖9 ωpet=100時系統(tǒng)電勢分布

3 結(jié) 論

在化學物質(zhì)釋放粒子模擬中,針對粒子模擬仿真模型,本文研究了泊松方程的FFT方法求解問題,提出一種改進的FFT方法求解泊松方程,取消近似替代項,并引入誤差修正項,消除非線性誤差對計算結(jié)果的影響,在仿真案例中求得的點電荷泊松方程電勢的MSE相對Birdsall方法和電勢求導方法小6個量級,誤差對數(shù)量級也遠遠小于這兩種方法,同時電勢空間分布更加符合實際．通過實際應用仿真模擬發(fā)現(xiàn),仿真結(jié)果符合理論推導,結(jié)果證明五點FFT方法求得電勢的精度符合化學物質(zhì)釋放早期小尺度不均勻體演化模擬的需要．

研究證明:改進的五點傅里葉變換方法在傅里葉空間的系數(shù)運算中并未使用近似替代等處理方法,減小了非線性誤差,同時對二階差分誤差采取了誤差彌補措施,進一步提高了差分精度,保證了算法的穩(wěn)定性,這種算法精度高、穩(wěn)定性好,而且計算效率尚可;采用改進的傅里葉變換方法求得電勢的空間分布更接近解析法求解結(jié)果,精度更高,能夠滿足化學物質(zhì)釋放數(shù)值仿真的精度要求,具有重要的理論意義和實際應用價值．在未來的工作中該方法將應用于中性氣體釋放產(chǎn)生電子耗空的混合模型數(shù)值模擬,以及電子增強的化學物質(zhì)釋放數(shù)值仿真中．

[1] TREFETHEN L N. Spectral methods in Matlab[M]. Philadelphia:Society for Industrial &Applied, 2000.

[2] 馬征征, 徐彬, 許正文, 等. 塵埃等離子體混合模型中電荷-電勢問題的迭代法求解[J]. 電波科學學報, 2015, 30(3): 549-553.

MA Z Z, XU B, XU Z W, et al. Solving charge-potential issue in dusty plasma hybrid model by iteration method[J]. Chinese journal of radio science, 2015, 30(3): 549-553. (in Chinese)

[3] KRUER W L, DAWSON J M, ROSEN B, et al. The dipole expansion method for plasma simulation[J]. Journal of computational physics, 1973, 13(1): 114-129.

[4] 邵福球. 等離子體粒子模擬[M]. 北京: 科學出版社, 2002.

[5] COOLEY J W, TUKEY J W. An algorithm for the machine calculation of complex Fourier series[J]. Math comput. 1965.

[6] WEAVER H J. Applications of discrete and continuous Fourier analysis [M]. New York: Wiley, 1983.

[7] BIRDSALL C K, LANGDON A B. Plasma physics via computer simulation[M]. New York: Mcgraw-Hill Book Company, 1985.

[8] SCALES W A, BERNHARDT P A, GANGULI G, et al. Early time evolution of negative ion clouds and electron density depletions produced during electron attachment chemical release experiments[J]. Journal of geophysical research, 1994, 99(A1): 373-381.

[9] CHAE G S, SCALES W A, CHAIR C, et al. Numerical simulation of ion waves in dusty plasmas[D]. Virginia: Virginia Polytechnic Institute and State University, 2000.

[10] 盧文祥.信號分析[M].武漢: 華中科技大學出版社.2002.

[11] 蔣伯誠.計算物理中的譜方法[M]. 湖南科學技術(shù)出版社, 1988.

[12] 馮象初, 任春麗,尚曉清等.數(shù)值分析[M]. 西安電子科技大學出版社.2012.