塵埃等離子體混合模型中電荷-電勢問題的迭代法求解

馬征征 徐 彬 許正文

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

引 言

在電離層、等離子體和電波傳播等領域的工程實際問題中,往往都需要求解電勢場分布. 對于靜電場,在電荷分布已知的條件下,通常由泊松方程來求解電勢. 由于泊松方程是偏微分方程,難以獲得解析解,因此，通常都使用數值計算方法. 常用的方法有解矩陣法、譜方法和迭代法等[1].

然而,許多空間等離子體環境中有可能出現塵埃粒子. 這些塵埃將不可避免地與背景等離子體相互作用:一方面,浸沒在等離子體中的塵埃粒子將通過碰撞和電離等方式充放電荷,成為帶電塵埃粒子;另一方面,帶電塵埃將顯著影響到等離子體中的動力學、電學、波動和不穩定等過程. 最終兩者相互作用,形成所謂塵埃等離子體. 常見的塵埃等離子體現象有地球極區夏季冰晶云、行星際空間、彗星、行星環以及火箭噴焰等[2-4].

帶電塵埃粒子的出現會使電勢場求解變得復雜. 首先,對于總電荷分布,除了電子和離子還需增加考慮塵埃粒子. 由于塵埃攜帶電荷量可以很大且本身具有很大的慣性,容易形成非常不均勻的電荷空間分布. 另由于塵埃粒子相比電子具有大得多的質荷比,在模擬和求解某些問題時,需要對兩者采用不同的處理方式,建立所謂混合模型. 常見的一類處理方式是:對塵埃粒子和離子采用動力學模擬,如質點網格(Particle-in-Cell,PIC)法,通過運動方程計算其分布;而對電子采用流模擬,如直接采用波爾茲曼分布[5-6]. 很重要地,此時電子電荷分布按照波爾茲曼分布律由電勢決定;而電勢通過泊松方程又由總電荷分布決定. 該問題不是直接由電荷求解電勢的過程,我們將其稱為“電荷-電勢”問題. 可以想到,對于此問題,解矩陣法和譜方法的適用性降低.

在如分子生物物理學等其他領域中,也存在著類似的問題,并且已開展了不少針對該問題中泊松-波爾茲曼方程求解方法的研究[7]. 而在塵埃等離子體領域,針對此問題的研究較少. 本文針對塵埃等離子體領域的工程實際,考察迭代法的求解適用性.

1 求解原理

首先給出求解靜電場中電勢的泊松方程為

(1)

式中:ρ是電荷密度; 下標d、i、e分別代表塵埃粒子、離子和電子;ε0是真空中的介電常數.

混合模型中,塵埃粒子電荷密度ρd和離子電荷密度ρi由動力學模型計算得到. 由于本文不研究動力學模型,因此直接作為已知量給出. 電子電荷密度ρe滿足波爾茲曼分布,有

(2)

式中:e是電子電量;kB是波爾茲曼常數;Te是電子溫度;ρe0是背景電子電荷密度.

最終由式(1)和(2)得到方程

f(φ) =2φ

=0.

(3)

注意到式(3)是偏微分的超越方程,解矩陣法和譜方法的適用性較低,因此考慮迭代法. 牛頓迭代法是平方收斂的,具有收斂速度快的特點. 其公式為

(4)

n表示迭代次數.

式(3)中還可以增加所謂下山因子ω,也稱為牛頓下山法. 格式為

(5)

若采用固定值A0替代f′(φn),則

(6)

則被稱為平行弦法,A0有時可取為f′(φ0). 下標0表示初始猜測值. 平行弦法的收斂速度通常相當于或低于牛頓迭代法,但省去了對函數導數的計算,節約了計算時間,且能夠比牛頓迭代法更加穩定. 其中,下山因子的確定有理論依據,通常以滿足|f(φn+1)|<|f(φn)|為原則. 但考慮到工程實際問題的復雜性,往往直接采用嘗試取值的方法. 此外下山因子還可隨迭代過程變化而取值[8].

2 求解過程及結果

2.1 牛頓下山法求解的簡單例

給出,且最小值為4.8×10-8C/m3. 其中x和y分別是網格橫向和縱向編號. 如圖1所示.

圖1 塵埃粒子和離子總電荷密度二維分布簡單例

電勢的初始值按均勻分布給出,通過式(5)來迭代求解. 迭代過程中期望f(φ)的值逐漸趨近于0,f(φ)=0對應的電勢值即為方程的解. 雖然實際求解過程中f(φ)=0往往難以實現,但f(φ)的值越接近0,對應的解也就越精確. 可用函數f(φ)的值來表征解的精度. 注意到f(φ)實際上是一個矩陣,可以通過定義式

(7)

來表征迭代過程中解的精度. 對于下山因子,考慮了四種不同值. 最終得到了迭代過程中f值(對數)的變化,如圖2所示.

圖2 不同下山因子下牛頓下山法求解簡單例的收斂過程

可以看到當下山因子ω=0.2、0.5和1.0時,迭代過程均收斂. 同時,收斂速度隨ω增大而加快. 當ω=1.0時,收斂速度非常快(9步). 此外,收斂后最終的f值約為10-13,表示解的精度非常高. 而當ω增大到2.0時,迭代不具收斂性. 這里我們取ω=1.0,得到迭代次數n=50時的電勢解如圖3所示,再由式(2)得到電子電荷密度解,如圖4所示.

2.2 平行弦法求解工程實際中的復雜例

考察來自工程實際中塵埃等離子體離子聲波不穩定性過程模擬的一個例子.網格為64×64,邊長為0.005 m,采用周期邊界條件[6]. 離子和帶電塵埃粒子的分布為已知量,已由動力學模型計算得到,如圖5所示. 可以看到,在帶電塵埃粒子加入時,電荷的空間分布可以十分復雜.

圖3 牛頓下山法求解簡單例的電勢二維分布

圖4 牛頓下山法求解簡單例的電子電荷密度二維分布

圖5 塵埃粒子和離子總電荷密度二維分布復雜例

圖6 不同下山因子下牛頓下山法求解復雜例的收斂過程

仍嘗試運用牛頓下山法求解,得到不同下山因子下f值隨迭代過程的變化結果,如圖6所示.可以看到當ω=0.01、0.1和1.0時,迭代過程在最初均表現出了收斂性,且收斂速度隨ω增大而加快. 然而,若干步后迭代均終止,表明牛頓下山法失敗. 分析表明,迭代若干步后f′(φn)矩陣中少數元素計算值過小,導致由式(5)計算的電勢值變化過大,直接違背了物理意義.從迭代過程的最初表現來看趨勢還是不錯的,因此，解決分母上的小值問題是關鍵. 可以考慮用固定值替代f′(φn),即所謂平行弦法. 可采用f′(φ0)來替代,有

(8)

利用平行弦法計算得到不同下山因子下f值的變化過程,結果如圖7所示.

圖7 不同下山因子下平行弦法求解復雜例的收斂過程

從圖7分析可知,下山因子在一定范圍內時(ω=0.1～0.5),平行弦法具有了收斂性,且迭代收斂速度隨ω值增大而加快. 但下山因子取值過大時(ω=0.8),求解出現問題. 圖8是ω=0.5、n=1 000時的電勢二維分布計算結果. 再由式(2)得到電子電荷密度二維分布,如圖9所示.

對于該復雜例,平行弦法的收斂速度是比較慢的. 圖7中的結果表明,迭代次數達到1 000步時,收斂仍未完成. 同時,從解精度的預期來看,也不如簡單例. 這是由兩例中離子和塵埃粒子總電荷密度空間分布的復雜性差異造成的.

圖8 平行弦法求解復雜例的電勢二維分布

圖9 平行弦法求解復雜例的電子電荷密度二維分布

3 結 論

塵埃等離子體領域工程實際問題中,經常采用混合模型處理不同的組分. 我們針對其中的一類情形研究了電荷-電勢問題的求解. 具體情形為電子密度分布由波爾茲曼分布律決定,而電勢場通過泊松方程由總電荷分布決定.

研究結果表明:對于簡單的情況,牛頓下山法完全適用,且具有非常快的收斂速度和非常高的求解精度;而對于復雜的工程實際例子,牛頓下山法不再適用. 分析研究表明,平行弦法對此類工程實際問題具備可觀的適用性. 下山因子取值在一定范圍內時,迭代具備收斂性和穩定性,且收斂速度隨下山因子的增大而加快. 當然,進一步的研究工作將在后續開展.

[1] TREFETHEN L N. Spectral Methods in MATLAB [M]. Philadelphia: Society for Industrial & Applied Mathematics, 2000.

[2] SHUKLA P K, MAMUN A A. Introduction to Dusty Plasma Physics[M]. London: IOP Publishing Ltd, 2002.

[3] 陳惠敏, 石雁祥, 蘇金善, 等. 基于動力論的塵埃等離子體充電方程研究[J]. 電波科學學報, 2010, 25(5): 1005-1009.

CHEN Huimin, SHI Yanxiang, SU Jinshan, et al. Charging equation of dust plasmas based on kinetic theory[J]. Chinese Journal of Radio Science, 2010, 25(5): 1005-1009. (in Chinese)

[4] 石雁祥, 王 菊, 吳 健, 等. 對兩種弱電離塵埃等離子體特征參量的定量估計[J]. 電波科學學報, 2008, 23(1): 95-99.

SHI Yanxiang, WANG Ju, WU Jian, et al. Characteristic parameters estimation of two weakly ionized dusty plasma[J]. Chinese Journal of Radio Science, 2008, 23(1): 95-99. (in Chinese)

[5] WINSKE D, YIN L, OMIDI N, et al. Hybrid Simulation Codes: Past, Present and Future—A Tutorial[M]. Berlin: Springer-Verlag, 2003: 136-165.

[6] CHAE G S, SCALES W A, GANGULI G, et al. Numerical model for low-frequency ion waves in magnetized dusty plasmas[J]. IEEE Trans Plasma Sci, 2000, 28(5): 1694-1705.

[7] HOLST M J. The Poisson-Boltzmann Equation-Analysis and Multilevel Numerical Solution (Monograph based on the Ph.D. Thesis: Multilevel Methods for the Poisson-Boltzmann Equation)[R]. Los Angeles: California Institute of Technology, 1994.

[8] 鄭慧嬈, 陳紹林, 莫忠息, 等. 數值計算方法[M]. 武漢:武漢大學出版社, 2002.