數(shù)值差分算法在靜電場模擬實(shí)驗(yàn)中的設(shè)計(jì)與應(yīng)用

翟立朋，常凱歌*，程 琳，張俊武，童 童

(1.西安交通大學(xué) 物理學(xué)院 物理國家級(jí)實(shí)驗(yàn)教學(xué)示范中心，陜西 西安 710049；2.西安工程大學(xué) 理學(xué)院，陜西 西安 710048)

模擬法描繪靜電場是普通物理實(shí)驗(yàn)中電磁學(xué)模塊的一個(gè)核心實(shí)驗(yàn)。帶電體在其周圍產(chǎn)生的靜電場可以通過電場強(qiáng)度E或電勢φ的空間分布來描述。對(duì)類似于同軸圓柱、圓筒這樣的規(guī)則電極，可以根據(jù)已知的電荷分布用靜電場方程求得電勢和電場分布的解析解。而對(duì)于不規(guī)則電極來說，電勢和電場分布的數(shù)學(xué)求解十分困難，需要用實(shí)驗(yàn)方法進(jìn)行測量[1]。

模擬法本質(zhì)上是通過一種容易實(shí)現(xiàn)、便于測量的物理量和物理過程，去代替不容易實(shí)現(xiàn)的、不便于測量的物理量和物理過程[2,3]。靜電場模擬實(shí)驗(yàn)通過測量相對(duì)容易的恒定電流去實(shí)現(xiàn)對(duì)靜電場的測量，因?yàn)殪o電場中沒有電流不能使用電磁式儀表進(jìn)行直接測量，同時(shí)探測裝置必須是導(dǎo)體或電介質(zhì)，一旦將其放入靜電場中，將會(huì)產(chǎn)生感應(yīng)電荷引起電場畸變，影響測量結(jié)果的準(zhǔn)確性。恒定電流之所以能模擬靜電場是因?yàn)殡娏鱆和電場E有相同的物理表達(dá)形式[4]。近年來，針對(duì)模擬法測量靜電場實(shí)驗(yàn)的數(shù)據(jù)處理方法[5,6]、儀器改進(jìn)[7-9]、誤差分析[10,11]、教學(xué)方法[12，13]與測量手段[14,15]的創(chuàng)新等都已經(jīng)有了很好的研究，但實(shí)驗(yàn)?zāi)Ｐ腿匀痪窒抻谘芯客S圓柱、圓筒電極間的靜電場。在普通物理實(shí)驗(yàn)上對(duì)于不規(guī)則電極，一般只測量其等勢線，再根據(jù)等勢線與電場線的關(guān)系得到電場線的分布定性的畫出電場線，無法真正得到靜電場的具體數(shù)值。

隨著計(jì)算機(jī)信息技術(shù)的發(fā)展，數(shù)值方法越來越多地被應(yīng)用到教學(xué)和科研中來。實(shí)驗(yàn)也可以使用COMSOL、Ansys、FDTD Solution等成熟的商業(yè)軟件進(jìn)行仿真分析[16,17]，但商業(yè)軟件往往費(fèi)用昂貴并且學(xué)生在使用商業(yè)軟件的過程中往往無法完整理解其仿真算法，所以不太適合在大面積基礎(chǔ)物理實(shí)驗(yàn)課程中進(jìn)行推廣。基于MATLAB將數(shù)值差分算法應(yīng)用到靜電場模擬實(shí)驗(yàn)中，可以給出任意形狀電極靜電場、電荷的分布圖，一方面擴(kuò)展了該實(shí)驗(yàn)項(xiàng)目的應(yīng)用場景，更重要的是將科研思維和方法應(yīng)用到物理實(shí)驗(yàn)教學(xué)中，可以大面積在基礎(chǔ)物理實(shí)驗(yàn)課程中推廣應(yīng)用。

1 模擬靜電場的基本原理

由電磁場理論可知[4]，均勻?qū)щ娊橘|(zhì)中穩(wěn)恒電流所產(chǎn)生的電流場與均勻電介質(zhì)中的靜電場具有相似性，這兩種物理場所遵循的物理規(guī)律具有相同的數(shù)學(xué)形式。對(duì)于靜電場而言，當(dāng)介質(zhì)中沒有自由移動(dòng)的電荷時(shí)，電場強(qiáng)度矢量E滿足方程：

(1)

∮lE·dl=0

(2)

對(duì)于穩(wěn)恒電流，當(dāng)介質(zhì)中無電流源時(shí)，電流場中的電流密度矢量J滿足方程：

(3)

∮lJ·dl=0

(4)

由此可見E和J在各自區(qū)域中滿足同樣的數(shù)學(xué)規(guī)律[1-3]。并且這兩個(gè)物理量都與電勢φ之間有著內(nèi)在的關(guān)聯(lián)，所以在實(shí)驗(yàn)上由直流電路構(gòu)建一個(gè)穩(wěn)恒的電流場，通過測量電極之間電勢差得到電流場的分布，進(jìn)而替代靜電場的分布。

圖1 同心圓柱、圓筒電極的實(shí)物與模型圖

如圖1所示，為無限長且同軸的圓柱和圓筒電極的橫截面，圓柱A和圓筒B分別作為正、負(fù)電極，并設(shè)定圓筒B的電勢為零，即φB=0。rA和rB分別表示圓柱和圓筒的半徑，如果以圓柱A的圓心為坐標(biāo)原點(diǎn)，則A、B之間位于r處的電勢可以通過高斯定理[3]得到：

(5)

靜電場是電勢的負(fù)梯度，電勢φr對(duì)位置r取導(dǎo)數(shù)，就可以得到任意位置處的電場強(qiáng)度表達(dá)式：

(6)

所以，只需要將A、B兩端電極加上穩(wěn)定的電壓，測量電路中不同位置處的電勢即可得到φr的具體數(shù)值，根據(jù)位置信息可以計(jì)算得到具體的電場強(qiáng)度值。

為了滿足上述理論中的需求，在實(shí)驗(yàn)上使用如圖2所示的電路在A、B兩個(gè)電極之間加一個(gè)可調(diào)電勢差。電極之間的灰色矩形框代表導(dǎo)電微晶，導(dǎo)電微晶是一種電阻分布均勻的導(dǎo)體，移動(dòng)電勢探頭的位置，電壓表的讀數(shù)即為電極B與探頭之間的電勢差，如果設(shè)定電極B的電勢為零，則電壓表可以測量任意位置處的電勢。

圖2 模擬靜電場實(shí)驗(yàn)電路圖

但是遺憾的是，實(shí)驗(yàn)雖然可以計(jì)算得到電場強(qiáng)度，但是得到電場強(qiáng)度的數(shù)值的方法與各個(gè)位置處的電勢沒有關(guān)系，在教學(xué)過程中體現(xiàn)不出模擬法具體應(yīng)用；更重要的是實(shí)驗(yàn)雖然可以測量任意形狀、任意位置處的電勢，但在理論上只能對(duì)同軸的圓柱、圓環(huán)模型進(jìn)行解析計(jì)算，其他形狀的電極無法給出類似于公式(6)這樣的解析表達(dá)式。所以，設(shè)計(jì)并應(yīng)用數(shù)值差分方法，將實(shí)驗(yàn)推廣到任意電極情況下的靜電場模擬。

2 數(shù)值差分算法及可行性驗(yàn)證

如果可以得到電勢隨位置的變化關(guān)系圖像，根據(jù)電磁學(xué)理論可以計(jì)算電場強(qiáng)度的大小。對(duì)于任何形狀的電極，電場強(qiáng)度為電勢的負(fù)梯度，這是恒定成立的，即：

E=-?φ

(7)

對(duì)于如公式(5)所示，電勢有解析表達(dá)式的電極模型，通過公式(7)可以計(jì)算得到電場強(qiáng)度的解析表達(dá)式。而對(duì)于不規(guī)則電極，引入數(shù)值差分算法對(duì)公式(7)進(jìn)行求解。

在上述同心圓柱、圓筒模型中，電勢沿著徑向成空間一維分布。對(duì)于梯度差分來說，數(shù)值差分算法的第一步是將位置空間進(jìn)行離散化，所以首先需要的是將位置r每間隔Δr進(jìn)行劃分，這樣空間的任意位置都可表示為：

r=n·Δr

(8)

其中n非負(fù)整數(shù)，根據(jù)對(duì)空間位置的離散劃分，任意位置處的電勢可表示為：

φr=φ(n·Δr)

(9)

在完成第一步空間離散化以后，第二步需要將電勢對(duì)空間的微分近似為差分，電勢的梯度表示為：

(10)

微分與差分的區(qū)別在于極限Δr的取值精度，Δr取值越小極限精度越高。所以，只要能保證離散空間Δr足夠小，就可以保證差分算法的計(jì)算結(jié)果足夠精確。當(dāng)然Δr不能無窮小，否則也無法進(jìn)行計(jì)算機(jī)的仿真計(jì)算?？梢圆痪窒抻谔囟ǖ碾姌O模型，將任意位置處的電場強(qiáng)度表示為：

(11)

這樣不需要解析解也可以描述出任意位置處的電場強(qiáng)度。

為了驗(yàn)證該方法的可行性，使用GVZ-4型靜電場描繪儀對(duì)同軸圓柱、圓筒模型進(jìn)行測量和計(jì)算，分別通過解析算法和數(shù)值差分算法對(duì)電勢、電場的測量和計(jì)算結(jié)果進(jìn)行對(duì)比分析。調(diào)節(jié)電路使電極A、B的電勢分別為φa=10.00 V、φB=0.00 V。該模型中同心圓柱的半徑為ra=1.0 cm，同心圓環(huán)的內(nèi)徑為rB=6.4 cm，將參數(shù)帶入到公式(5)和公式(6)可以計(jì)算得到A、B之間任意位置處電勢的解析表達(dá)式為：

φr=10.00-5.39·lnr

(12)

可以看出φr與lnr之間滿足線性關(guān)系。得到A、B之間任意位置處電場強(qiáng)度的解析表達(dá)式為：

(13)

實(shí)驗(yàn)上用圖1所示的電極模型，測量過程中使用的位置差分間隔為Δr=0.2 cm，將測量得到的電勢值代入公式(11)得到對(duì)應(yīng)位置處的電場強(qiáng)度值，測量得到各個(gè)半徑r處對(duì)應(yīng)的電勢和電場強(qiáng)度數(shù)值如表1所示。需要說明的是電場強(qiáng)度是相鄰兩個(gè)電勢做差之后計(jì)算得到，差分之后的電場強(qiáng)度數(shù)據(jù)量要比原始的電勢數(shù)據(jù)少一個(gè)，所以表1中r=1.2 cm處沒有電場強(qiáng)度值。

表1 同軸圓柱圓筒電勢與電場值

如圖3(a)和(b)分別表示同軸圓柱圓筒模型中電勢和電場強(qiáng)度值隨著位置的變化關(guān)系，圖像中的直線與曲線代表的是理論值，圖3(a)中的散點(diǎn)代表的是電勢的測量值，圖3(b)中的散點(diǎn)代表的是經(jīng)過差分公式計(jì)算得到的電場值。

lnr/cm(a)

可以看出電勢測量值以及經(jīng)過差分計(jì)算的電場值與理論值吻合非常好。

進(jìn)一步計(jì)算得到電勢測量值的標(biāo)準(zhǔn)偏差為：

(14)

其中φi和φ0分別表示對(duì)應(yīng)位置處電勢的測量值與理論值。差分結(jié)果與理論值之間的標(biāo)準(zhǔn)偏差為：

(15)

圖3所展示出的關(guān)系圖以及標(biāo)準(zhǔn)偏差，都說明了該方法在模擬靜電場實(shí)驗(yàn)中是可行的。

對(duì)于其他形狀的電極來說，需要建立二維空間直角坐標(biāo)系，可以通過如圖2所示的電路測量得到二維空間的電勢分布φ(x,y)。二維空間離散化以后任意空間位置可以表示為x=m·Δx、y=n·Δy，則電場強(qiáng)度的分布就可以表示為：

(16)

(17)

在得到電場分量得到以后即可得到電場強(qiáng)度的絕對(duì)值：

(18)

該算法不僅可以計(jì)算任意形狀的電場分布，進(jìn)行二次差分以后，還可以得到空間電荷分布。根據(jù)高斯定理，空間任意位置處的電荷分布可表示為：

(19)

(20)

(21)

將(20)、(21)兩式代入到公式(19)以后，就可以得到空間任意位置處的電荷分布，這樣只要測量得到任意電極模型的電勢分布，靜電場的電場強(qiáng)度、電荷等參量即可全部得到。

3 任意形狀電極電場與電荷分布

如圖4所示，為GVZ-4型靜電場描繪儀中的一種四電極模型，并建立如圖所示的xoy坐標(biāo)系，坐標(biāo)系中導(dǎo)電微晶上每一個(gè)最小網(wǎng)格為：Δx=0.2cm、Δy=0.2cm，以此為離散空間的大小。下面兩個(gè)電極都接電路中的B電極，這兩個(gè)電極的電勢都為φB=0.00V，上面兩個(gè)電極都為A電極并設(shè)定電勢為φa=10.00V。

圖4 四電極模型的結(jié)構(gòu)圖與坐標(biāo)系

每間隔0.2cm測量各個(gè)位置上的電勢以后，就可以得到如圖5所示的電勢分布圖，圖中每條曲線及對(duì)應(yīng)的數(shù)值代表等勢線及對(duì)應(yīng)的電勢值。從等勢線可以清晰地看出電極輪廓及二維平面上的電勢分布。因?yàn)殡姌O都是金屬導(dǎo)體，從等勢線也可以看出四個(gè)電極都為等勢體，下面兩個(gè)電極電勢為0.00V，上面兩個(gè)電極電勢為10.00V，這都符合實(shí)驗(yàn)設(shè)定的初始條件。

圖5 四電極模式的電勢分布

測量得到二維電勢分布以后，使用MATLAB對(duì)二維數(shù)據(jù)進(jìn)行差分運(yùn)算，MATLAB差分計(jì)算的整體模型和流程如圖6所示。

圖6 MATLAB數(shù)值差分計(jì)算模型

分別將二維電勢分布在x方向和y方向進(jìn)行差分和插值，這里需要插值補(bǔ)齊是因?yàn)榻?jīng)過差分運(yùn)算以后會(huì)有數(shù)據(jù)損失，例如在表1中電勢值有27個(gè)數(shù)據(jù)，但是差分之后r=1.2cm處就少了一個(gè)電場數(shù)據(jù)，電場強(qiáng)度只有26個(gè)數(shù)據(jù)，需要經(jīng)過插值運(yùn)算將26個(gè)數(shù)據(jù)補(bǔ)齊到27個(gè)，計(jì)算得到的電場和電荷分布都使用的是三次樣條插值算法。

如圖7所示的電場分布可以看出，強(qiáng)電場主要分布在A、B電極之間，同A電極或同B電極之間的電場強(qiáng)度較小。電極附近尖銳位置處的電場強(qiáng)度較強(qiáng)，符合尖端效應(yīng)的預(yù)期。

圖7 四電極模型電場強(qiáng)度分布

圖8所示的電荷分布中紅藍(lán)色分別代表正負(fù)電荷，白色區(qū)域代表無電荷分布，可以看出大量電荷主要分布在電極導(dǎo)體表面以及導(dǎo)體尖端處，導(dǎo)體內(nèi)部無電荷分布與靜電屏蔽原理吻合，導(dǎo)電微晶平面上分布有少量負(fù)電荷，每兩個(gè)同A電極或同B電極之間的正負(fù)電荷都成對(duì)稱分布，說明這是在同電極之間兩個(gè)不同位置處形成了偶極子。

圖8 四電極模型電荷分布

4 結(jié) 語

將數(shù)值差分算法引入到模擬靜電場實(shí)驗(yàn)中，實(shí)現(xiàn)了電場強(qiáng)度分布和電荷分布的測量，彌補(bǔ)了該實(shí)驗(yàn)項(xiàng)目沒有完全測到電場數(shù)值的欠缺，雖然只給出了一種四電極模型的電場、電荷分布情況，但方法可以應(yīng)用到任意形狀的電極模型，這對(duì)實(shí)驗(yàn)內(nèi)容的深度和廣度是非常有效地拓展。