馬征征 徐 彬 許正文

(中國(guó)電波傳播研究所 電波環(huán)境特性及模化技術(shù)重點(diǎn)實(shí)驗(yàn)室,山東 青島 266107)

引 言

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

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

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

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

1 求解原理

首先給出求解靜電場(chǎng)中電勢(shì)的泊松方程為

(1)

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

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

(2)

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

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

f(φ) =2φ

=0.

(3)

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

(4)

n表示迭代次數(shù).

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

(5)

若采用固定值A(chǔ)0替代f′(φn),則

(6)

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

2 求解過(guò)程及結(jié)果

2.1 牛頓下山法求解的簡(jiǎn)單例

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

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

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

(7)

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

圖2 不同下山因子下牛頓下山法求解簡(jiǎn)單例的收斂過(guò)程

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

2.2 平行弦法求解工程實(shí)際中的復(fù)雜例

考察來(lái)自工程實(shí)際中塵埃等離子體離子聲波不穩(wěn)定性過(guò)程模擬的一個(gè)例子.網(wǎng)格為64×64,邊長(zhǎng)為0.005 m,采用周期邊界條件[6]. 離子和帶電塵埃粒子的分布為已知量,已由動(dòng)力學(xué)模型計(jì)算得到,如圖5所示. 可以看到,在帶電塵埃粒子加入時(shí),電荷的空間分布可以十分復(fù)雜.

圖3 牛頓下山法求解簡(jiǎn)單例的電勢(shì)二維分布

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

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

圖6 不同下山因子下牛頓下山法求解復(fù)雜例的收斂過(guò)程

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

(8)

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

圖7 不同下山因子下平行弦法求解復(fù)雜例的收斂過(guò)程

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

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

圖8 平行弦法求解復(fù)雜例的電勢(shì)二維分布

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

3 結(jié) 論

塵埃等離子體領(lǐng)域工程實(shí)際問(wèn)題中,經(jīng)常采用混合模型處理不同的組分. 我們針對(duì)其中的一類情形研究了電荷-電勢(shì)問(wèn)題的求解. 具體情形為電子密度分布由波爾茲曼分布律決定,而電勢(shì)場(chǎng)通過(guò)泊松方程由總電荷分布決定.

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

[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] 陳惠敏, 石雁祥, 蘇金善, 等. 基于動(dòng)力論的塵埃等離子體充電方程研究[J]. 電波科學(xué)學(xué)報(bào), 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] 石雁祥, 王 菊, 吳 健, 等. 對(duì)兩種弱電離塵埃等離子體特征參量的定量估計(jì)[J]. 電波科學(xué)學(xué)報(bào), 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] 鄭慧嬈, 陳紹林, 莫忠息, 等. 數(shù)值計(jì)算方法[M]. 武漢:武漢大學(xué)出版社, 2002.