孫毅祥

(1. 北京長(zhǎng)征天民高科技有限公司， 北京 100176； 2. 中國(guó)運(yùn)載火箭技術(shù)研究院第十五研究所， 北京 100076)

0 引言

靜電學(xué)在工程中有諸多應(yīng)用， 如靜電場(chǎng)探測(cè)系統(tǒng)[1]、靜電屏蔽技術(shù)[2]等. 在靜電場(chǎng)的計(jì)算模擬中， 電勢(shì)求解是關(guān)鍵的一環(huán). 由于均勻電介質(zhì)中靜電場(chǎng)高斯定理可以轉(zhuǎn)換為泊松方程[3]， 因此靜電場(chǎng)不用求解復(fù)雜的麥克斯韋方程， 而只需求解二階橢圓型偏微分方程. 泊松方程常見(jiàn)的數(shù)值求解方法包括有限差分法[4]、譜方法[5]、有限元法等.

譜方法相比于有限差分法、有限元法而言， 其特點(diǎn)之一是用于近似的基函數(shù)不是局部點(diǎn)的低階多項(xiàng)式， 而是全局多項(xiàng)式[6](如Legendre 多項(xiàng)式、Chebyshev 多項(xiàng)式). 當(dāng)方程解光滑時(shí)其收斂速度很快， 可以在較少的網(wǎng)格點(diǎn)上實(shí)現(xiàn)高精度數(shù)值模擬[7，8]， 適合于泊松方程的高效準(zhǔn)確求解. 譜方法一般分為三類[9]： 譜-Galerkin法、譜配置法和譜-τ 法. 其中譜配置法將偏微分方程轉(zhuǎn)換為矩陣方程， 且對(duì)于Chebyshev 多項(xiàng)式其配置點(diǎn)可顯式給出， 因此無(wú)論從計(jì)算效率還是編程復(fù)雜度而言都是較優(yōu)的選擇. 對(duì)于泊松方程來(lái)說(shuō)， 若在極坐標(biāo)系下求解， 譜方法可以直接給出電勢(shì)在徑向和周向的變化， 但需要額外處理方程在極坐標(biāo)系原點(diǎn)奇異的問(wèn)題. 目前主要有兩種解決方法， 一種是構(gòu)造奇點(diǎn)處的特殊邊界條件[10，11]， 另一種則是選取合理的配置點(diǎn)(如Chebyshev-Gauss-Radau 配置點(diǎn))來(lái)避開(kāi)奇點(diǎn)位置[12，13]. 本文采用第二種奇點(diǎn)處理方式， 建立了極坐標(biāo)系下的泊松方程譜方法求解器.

雖然對(duì)均勻電介質(zhì)而言靜電場(chǎng)高斯定理可以轉(zhuǎn)化為泊松方程， 但對(duì)于非均勻電介質(zhì)， 原方程將多出介電系數(shù)的空間梯度項(xiàng)， 若此時(shí)仍用原始泊松方程進(jìn)行求解， 將會(huì)引入誤差， 且誤差隨介電系數(shù)的空間梯度增大而增大. 非均勻電介質(zhì)在實(shí)際中也有對(duì)應(yīng)情況， 如存在介質(zhì)分界面(水-空氣)或介質(zhì)本身連續(xù)非均勻. 因此電介質(zhì)非均勻的情形無(wú)論是理論研究還是實(shí)際應(yīng)用都有一定價(jià)值. 然而目前較少有研究直接關(guān)注非均勻介電系數(shù)對(duì)電勢(shì)的具體影響機(jī)制， 盡管近年有部分研究者注意到對(duì)修正泊松方程進(jìn)行求解更符合實(shí)際(如Nagel[14]的工作)， 但并未對(duì)介質(zhì)非均勻如何影響電勢(shì)分布進(jìn)行研究與分析. 本文基于譜方法求解器對(duì)非均勻電介質(zhì)的修正泊松方程進(jìn)行求解， 主要研究電介質(zhì)非均勻時(shí)仍采用原始泊松方程求解而引入的結(jié)果誤差， 以及電介質(zhì)非均勻?qū)Ψ匠探獾挠绊?

1 極坐標(biāo)系下修正泊松方程的求解

靜電場(chǎng)高斯定理的微分形式為

其中φ 為電勢(shì)， ρ 為電荷密度， ε 為介電系數(shù). 若在真空中或電介質(zhì)空間分布均勻， 介電系數(shù)可從梯度算子中提出， 式(1)便化為標(biāo)準(zhǔn)的泊松方程. 但對(duì)于電介質(zhì)非均勻的情況， 式(1)便轉(zhuǎn)化為

式(2)等號(hào)左端第二項(xiàng)即為電介質(zhì)非均勻時(shí)的額外項(xiàng)， 若忽略此項(xiàng)并按原始泊松方程進(jìn)行求解， 將會(huì)引入誤差. 下面推導(dǎo)修正泊松方程(2)在極坐標(biāo)系下的形式， 并給出譜配置法的求解細(xì)節(jié).

極坐標(biāo)系下， 梯度算子可寫作徑向和周向單位矢量的組合：

將其代入式(2)并在等號(hào)兩側(cè)乘以r2， 于是式(2)轉(zhuǎn)化為

其中包含g1、g2的項(xiàng)分別為介電系數(shù)在徑向和周向梯度引入的額外項(xiàng)， 表達(dá)式為

為方便使用譜配置法， 本文只考慮g1與g2分別僅與r 和θ 有關(guān)的情況， 即 g1= g1( r)， g2= g2(θ ). 此時(shí)式(4)左端系數(shù)已實(shí)現(xiàn)變量分離， 可將解φ 展開(kāi)為基函數(shù)與各配置點(diǎn)處解的乘積的線性組合. 本文考慮在去除極點(diǎn)的圓盤 (0， R0] × [0，2π )上求解式(4)， 徑向采用Chebyshev 多項(xiàng)式展開(kāi)， 同時(shí)為避開(kāi)奇點(diǎn)， 采用Chebyshev-Gauss-Radau 配置點(diǎn)

在周向則采用Fourier 展開(kāi)， 配置點(diǎn)為

注意到r 與配置點(diǎn) jα 定義域不同， 因此采用坐標(biāo)變換

此時(shí)式(4)轉(zhuǎn)化為

在徑向與周向分別展開(kāi)， 式(9)最終表示為矩陣方程的形式[15]：

系數(shù)矩陣及源項(xiàng)矩陣的元素為

其中徑向一階導(dǎo)數(shù)算子對(duì)應(yīng)的離散矩陣元素

這里

徑向二階導(dǎo)數(shù)算子對(duì)應(yīng)的離散矩陣

周向一階、二階導(dǎo)數(shù)算子對(duì)應(yīng)的離散矩陣則跟Nθ的奇偶性有關(guān)， 這里給出Nθ為偶數(shù)時(shí)的離散矩陣， 奇數(shù)情況參見(jiàn)文[9]：

以上即為修正泊松方程的離散求解形式. 若要考慮 r = R0處的第一類邊界條件， 只需在式(10)兩側(cè)減去 r = R0處的源項(xiàng)值F(BC)即可. 在本文中， F(BC)的各元素可寫作[15]

對(duì)形如式(10)的矩陣方程， 可采用文[15]中的兩步算法進(jìn)行求解， 最終得到電勢(shì)φ 的離散解.

2 原始泊松方程對(duì)非均勻電介質(zhì)的電勢(shì)求解誤差研究

若采用原始泊松方程對(duì)非均勻電介質(zhì)的電勢(shì)進(jìn)行計(jì)算， 即忽略式(2)左端第二項(xiàng)， 則會(huì)引入誤差， 本節(jié)擬對(duì)此誤差進(jìn)行評(píng)估. 首先給出計(jì)算采用的參數(shù)條件. 為與解析解進(jìn)行比較， 對(duì)解 φ = exp( - r2)sinθ(單位： V)對(duì)應(yīng)的修正泊松方程進(jìn)行數(shù)值求解， 源項(xiàng)可通過(guò)將解析的電勢(shì)解代入式(2)求出. 這里需要給定介電系數(shù)的表達(dá)式以數(shù)值求出φ. 本文擬分別考慮介電系數(shù)在徑向和周向梯度下的影響， 并對(duì)徑向部分舍去高階項(xiàng)而只考慮線性相關(guān)項(xiàng)， 因此給定如下表達(dá)式：

其中真空中介電常數(shù) ε0= 8.854 × 10-12(F/m)， 正弦函數(shù)可使得介電系數(shù)ε 在周向?yàn)橹芷诤瘮?shù)， 常數(shù)C1和C2其中之一為0， 以滿足式(5)中g(shù)1與g2分別只與r 和θ 有關(guān)的條件. 取 C3= 3， 目的是在模擬的參數(shù)范圍內(nèi)保證 ε ＞ ε0以符合物理實(shí)際. 此外半徑 R0= 1m.

首先給出網(wǎng)格收斂性與程序準(zhǔn)確性的驗(yàn)證. 取 Nr= Nθ， C1= C2= 0， 并定義數(shù)值解與解析解之差的二范數(shù)為誤差：

圖1給出了誤差隨網(wǎng)格點(diǎn)數(shù) Nr+ Nθ的變化情況. 不難看出， 隨著網(wǎng)格點(diǎn)數(shù)的增多， 誤差呈指數(shù)級(jí)下降， 且僅需 Nr= Nθ=10時(shí)絕對(duì)誤差就已下降到 10-9量級(jí)， 即譜配置法僅需很少的網(wǎng)格點(diǎn)就可以達(dá)到足夠高的精度， 證明了此方法的優(yōu)越性， 也驗(yàn)證了本程序的準(zhǔn)確性. 另外， 誤差下降到一定程度后會(huì)出現(xiàn)反彈， 這是由計(jì)算機(jī)誤差累計(jì)引起的[16]， 但此時(shí)誤差已經(jīng)足夠小， 不會(huì)對(duì)最終結(jié)果造成明顯影響.

圖1 數(shù)值解與解析解的誤差隨網(wǎng)格點(diǎn)數(shù)的變化

隨后， 對(duì)原始泊松方程(式(2)左端移去第二項(xiàng))和修正泊松方程(式(2))分別求解， 方程右端源項(xiàng)由解析解代入修正泊松方程計(jì)算得到， 介電系數(shù)給定為式(17)， 離散點(diǎn)數(shù)取 Nr= 20， Nθ= 40. 圖2 給出了 C1= 10， C2= 0時(shí)的解析解、修正泊松方程解和原始泊松方程解. 可看出， 修正泊松方程解與解析解擬合良好， 但原始泊松方程解卻產(chǎn)生了偏離， 除 θ = 0°和180°外， 其余等勢(shì)線相對(duì)于解析解均向外擴(kuò)張， 且上下半平面分別越靠近θ = 90°和270°， 與解析解的差異越明顯， 而在靠近 θ = 0°和180°附近差異不明顯， 即呈現(xiàn)出各向異性的特性.

以上是介電系數(shù)僅在徑向分布不均的情況. 當(dāng)介電系數(shù)僅在周向分布不均時(shí)， 圖3 給出了解析解和C1=0， C2=0.1 時(shí)的原始泊松方程解(修正泊松方程解與解析解仍然吻合很好， 這里不再給出). 不難看出， 介電系數(shù)周向分布不均時(shí)， 用原始泊松方程得到的解的等勢(shì)線不僅發(fā)生了偏移， 還破壞了解析解關(guān)于Y軸反對(duì)稱的特點(diǎn). 此外， 為評(píng)估原始泊松方程解與解析解的誤差， 仍用式(18)進(jìn)行計(jì)算. 結(jié)果發(fā)現(xiàn) C1= 10， C2= 0和 C1= 0， C2= 0.1時(shí)， 原始泊松方程的誤差分別約為1.21 和19.07， 而修正泊松方程的誤差在 1 0-14量級(jí). 綜合可知， 介電系數(shù)在周向存在很小的梯度( C2= 0.1)， 就會(huì)使得用原始泊松方程求解造成比介電系數(shù)存在較大徑向梯度時(shí)( C1= 10)求解更大的誤差， 即介電系數(shù)徑向梯度影響更顯著.

圖2 C1=10， C2=0 時(shí)電勢(shì)的解析解(云圖)、 修正泊松方程解(實(shí)線)、原始泊松方程解(虛線)

圖3 C1=0， C2=0.1 時(shí)電勢(shì)的解析解(云圖)、 原始泊松方程解(虛線)

3 非均勻電介質(zhì)對(duì)電勢(shì)的影響

本部分采用修正泊松方程(式(2))對(duì)非均勻電介質(zhì)中的電勢(shì)分布進(jìn)行求解， 目的在于研究介電系數(shù)不均勻程度發(fā)生改變時(shí)電勢(shì)分布的變化， 從而表征電介質(zhì)非均勻效應(yīng)的影響. 與第2 部分不同的是， 本部分固定的是式(2)中的電荷密度ρ 的表達(dá)式， 該表達(dá)式由 φ = exp( - r2)sinθ， ε= C3ε0代入式(2)得到， 代表均勻電介質(zhì)下的電荷密度； 而變動(dòng)的則是ε 表達(dá)式(17)中的常數(shù)C1或C2， 從而解φ 也是變動(dòng)的， 于是可以評(píng)估C1或C2對(duì)解φ 的影響.

圖4 C2 = 0 時(shí)不同C1 取值對(duì)θ = 45°處電勢(shì)的影響及漸進(jìn)解

圖4 與圖5 分別給出了C2=0 時(shí)， 不同C1取值對(duì)θ = 45°和 θ = 90°處電勢(shì)沿徑向分布的影響. 不難看出， 介電系數(shù)徑向梯度對(duì)θ = 45°和θ = 90°處的影響規(guī)律是相似的， 當(dāng) C1＞ 0即介電系數(shù)沿徑向有正梯度時(shí)， 電勢(shì)下降， 且隨C1的增大而下降更多； 而當(dāng) C1＜ 0即介電系數(shù)沿徑向有負(fù)梯度時(shí)， 電勢(shì)減小. 此外， 當(dāng)C1趨于無(wú)窮時(shí)， 電勢(shì)分布有漸進(jìn)解， 這是因?yàn)楹瘮?shù) g1在C1趨于無(wú)窮時(shí)有極限， 代入式(2)可求出電勢(shì)的漸進(jìn)解.

圖6 給出了 C1= 0時(shí)， 不同C2取值對(duì)θ = 45°處電勢(shì)沿徑向分布的影響. 可看出， C2＞ 0時(shí)電勢(shì)升高， 反之則降低. 同時(shí)， 電勢(shì)的變化相對(duì)于 C1≠ 0， C2= 0時(shí)更明顯， 如 C1= 0， C2= 0.5時(shí)最大電勢(shì)可達(dá)到均勻介質(zhì)( C1= C2= 0)時(shí)的13 倍左右， 而 C1= 10， C2= 0時(shí)最大電勢(shì)與均勻介質(zhì)的情況差異很小. 這意味著周向的影響更為顯著. θ = 90°處電勢(shì)分布規(guī)律與圖6 相似， 這里不再給出. 此外， C2趨于無(wú)窮時(shí)， 電勢(shì)分布也有漸進(jìn)解， 但函數(shù) g2的極限在θ = 0°及180°時(shí)是奇異的， 因此數(shù)值計(jì)算未能給出此時(shí)漸進(jìn)解的分布.

圖6 C1=0 時(shí)不同C2 取值對(duì)θ = 45°處電勢(shì)的影響

4 結(jié)論

以往研究中靜電場(chǎng)電勢(shì)求解問(wèn)題主要是求解均勻電介質(zhì)的泊松方程. 本文針對(duì)非均勻電介質(zhì)的修正泊松方程進(jìn)行了譜方法求解并研究介電系數(shù)非均勻的影響. 首先建立了高效準(zhǔn)確的譜方法求解器， 在較少的網(wǎng)格點(diǎn)下數(shù)值解就能與理論解良好吻合. 隨后， 對(duì)采用原始泊松方程求解非均勻電介質(zhì)的電勢(shì)分布而引入的誤差進(jìn)行了研究， 發(fā)現(xiàn)介電系數(shù)存在徑向和周向梯度時(shí)均會(huì)產(chǎn)生誤差， 其中有周向梯度時(shí)仍用原始泊松方程求解不僅對(duì)電勢(shì)數(shù)值的影響更大， 還會(huì)破壞精確解本身的分布特性(如關(guān)于Y 軸的反對(duì)稱特性). 最后， 用修正泊松方程準(zhǔn)確求解了非均勻電介質(zhì)的電勢(shì)分布， 并研究了介電系數(shù)的梯度對(duì)電勢(shì)分布的影響. 結(jié)果表明， 介電系數(shù)梯度越大， 電勢(shì)分布變化更大， 且周向梯度對(duì)電勢(shì)分布的影響更為顯著. 此外， 當(dāng)介電系數(shù)的徑向梯度和周向梯度趨于無(wú)窮時(shí)， 電勢(shì)分布具有漸進(jìn)解. 這些結(jié)論可以為靜電場(chǎng)電介質(zhì)非均勻效應(yīng)的后續(xù)分析及建模提供參考與啟發(fā).