Python求解不同形狀之電容值

Python 求解不同形狀之電容值

吳致頡* 姜理元 賴奕帆

臺北市立建國高級中學 *通訊作者：jackwu9923@gmail.com （投稿日期：民國106 年 12 月 05 日，接受日期：107 年 05 月 06 日） 摘要：本文簡介了有限差分法，應用於求解電位的拉普拉斯方程，其目的在於推 廣此物理中簡易卻實用的程式模擬技巧。文章中先說明有限差分法的基本原理， 諸如泰勒展開近似、網格、三種邊界條件等，以平行板電容器、圓柱形電容器和 球形電容器三種邊界條件為例，比較理論值及模擬值，最後求解一非標準形狀的 電容值，展現此方法之所以實用之處。 關鍵詞：有限差分法、拉普拉斯方程、電容、邊界值問題

壹、 前言

物理的研究越來越仰賴程式的模擬，許多問題在做實驗之前，若能有效模擬現象，對於 了解現象及實驗變因如何設定等疑問將極有幫助。近年來，許多商用軟體如 Wolfram Mathematica 和 Matlab，能方便的算出一些微分方程的解，畫出精美的圖和動畫，但這些軟體

通常要花錢購買，且許多人並不知道其中運作的道理。本文以python 結合 numpy 和 matplotlib

兩個免費函式庫，前者用於科學計算，後者用於繪圖，以有限差分法模擬滿足拉普拉斯方程 ，給定邊界條件下，不同形狀下之電容值。

希望推廣簡易且應用面廣的數值分析方式，能夠對大家無論是做學問或教學上都能有所 幫助。

貳、 理論

一、 有限差分

電腦的世界中，必須把現實世界幾近連續的事物拆成離散的才能計算。而這個簡化的過 程，最簡單的便是「差分」。而有限差分法又分為顯示法（explicit method）和隱式法（implicit method），前者的概念較容易，也就是本文將介紹的，而後者則須將方程式中的係數寫為一個 矩陣，以反矩陣求解，理論較為複雜，在此不介紹。 (一) 微分 1. 一次微分 由泰勒展開式： y = f(x) = ∑𝑓(𝑘)(x) 𝑘! (𝑥 − 𝑎) ∞ 𝑘=0 f(𝑥 + ℎ) ≈ f(𝑥) + f′_(x)h 故 f′_{(x) ≈}f(𝑥 + ℎ) − f(𝑥) h 差分中習慣以下標表示其座標，以 y= f(𝑥)為例，上式即可被表示為： dy dx≈ 𝑦𝑖+1−𝑦𝑖 ∆x 此稱為一階精度（只取泰勒展開的第一項）之差分。在二維的情形中，則表示為∂y_∂x≈ 𝑦𝑖+1,𝑗−𝑦𝑖,𝑗 ∆x ，以此類推，其離散的方式可分為以下三種：

(1) 前方差分（forward difference method） dy dx≈ ∆y ∆x= 𝑦_𝑖+1−𝑦_𝑖 ∆x (2) 中央差分（central difference method）

dy dx≈ ∆y ∆x= 𝑦𝑖+1−𝑦𝑖−1 2∆x (3) 後方差分（backward difference method）

dy dx≈ ∆y ∆x= 𝑦𝑖−𝑦𝑖−1 ∆x 一般而言，中央差分具有較高的精準度，較常用於計算在空間中的場。而若系統中的 數值隨時間變化，即t 為函數的變數時，對 t 通常採前方差分，因我們並不知道欲計 算的這個時間點上，函數的值為何，技術上無法將這個時間點與下一個時間的的值納 入運算。

一些數值分析中，雖然不考慮系統隨著時間的變化，但在求解的過程中，會不斷以前 一個值來疊代（iteration），以逼近真實值。疊代越多次，結果就越準確，最後值會趨 於一個定值。這樣的方法稱為relaxation method，而一次又一次的疊代就好像數值隨 著時間演變一樣，故被稱虛擬時間(pseudo time)。 本文模擬皆使用空間中央差分、虛擬時間前方差分法。 2. 二次微分 用類似的想法，可以寫成： 𝜕2_y 𝜕𝑥2≈ ∆ (∆y_∆x) ∆x = (𝑦𝑖+1,𝑗_∆x−𝑦𝑖,𝑗) − (𝑦𝑖,𝑗−𝑦_∆x𝑖−1,𝑗) ∆x = 𝑦_𝑖+1,𝑗−2𝑦_𝑖,𝑗+ 𝑦_{𝑖−1,𝑗} ∆x2 (二) 網格 1. 基本原理 為了將上述的差分方法應用在程式中，我們得視情況設定「網格」，也就是建立一個 空間來儲存上述方法所算出來的資料，而矩形是最簡單網格的樣子。 圖1：本程式用的網格 本程式用的是每格長寬（dx 和 dy）皆相同的網格，Ex和Ey是由電位計算出電場之x 與y 分量，並不納入計算；在此畫在線與線之間，因其為網格兩點間「平均變化率」 的概念。 2. 積分 有了網格的概念，不只能運用於求解微分方程，也可以用最基礎黎曼和的概念，對網 格內的數值積分。舉前一點為例，若要求點場平方對空間的積分（假設總寬度為 W 、高度為L，此面有深度 B（m），且電位值不隨深度而變），則可以寫出： ∭ 𝐸2_{𝑑𝑥𝑑𝑦𝑑𝑧 ≈ 𝐵 ∑ ∑ (𝐸} 𝑥𝑖,𝑗2+ 𝐸𝑦𝑖,𝑗2) 𝑛𝑥−1 𝑖=0 𝑛𝑦−1 𝑗=0

如此便能得到積分的結果。在python 中，這樣的式子可以寫成： for i in range(nx):

for j in range(ny):

integration += sum((Ex[j,i]**2 + Ey[j,i]**2)*(dx*dy*B)) 或者，較為簡便的寫法：integration = (Ex**2 + Ey**2)*(dx*dy*B) (三) 邊界條件

微分方程須搭配邊界條件方能得解，分為三種，即

(1) 第一類邊界條件（即，狄利克雷邊界條件，Dirichlet boundary condition）

已知所求解之微分方程中，包圍此系統所有邊界的物理量。例如：一個圓柱形電 容中，已知內外圓柱的電位。欲求解的區域在兩圓之間，被完全包圍。

(2) 第二類邊界條件（即，諾伊曼邊界條件，Neumann boundary condition）

已知所求解之微分方程中所有邊界上，物理量對於邊界垂直方向的導數。程式中， 這樣的條件可以寫成離散的𝑦𝑖+1 = 𝑦𝑖+ 𝑦′𝑑𝑥。本文中未使用此種邊界條件，但是 系統如一座湖中央的一塊區域（二維），已知這塊區域四周垂直區域邊緣的水流流 速，便可由第二類邊界條件求解這一塊區域內各處的流速。 (3) 第三類邊界條件（混合型邊界條件） 已知所求解之微分方程中所有邊界上，知道物理量或物理量對於邊界垂直方向導 數的其中之一。例如：平行電容板（側視），上下板的電位已知，而左右方的缺口 雖然不知道電位為何，但已知左右邊界上電場（電位的導數）沒有水平分量。所 有已知的條件也把欲求解的區域完全包圍起來了。 對於拉普拉斯方程而言，當包圍此空間的邊界條件都已確定（無論是上述的哪一個） ，其解必唯一存在。

二、 電學

(一) 電位的 Laplacian=0 𝐸⃗ = −∇V 由高斯定律的微分形式： ∇ ∙ 𝐸⃗⃗⃗⃗ = 𝜌 𝜀₀ 因為空氣為電中性，ρ=0，故： ∇2_{V =}𝜕2V 𝜕2_x+ 𝜕2_V 𝜕2_y+ 𝜕2_V 𝜕2_z= 𝜌 𝜀0= 0 --式 1 或圓柱座標的形式：

∇2_{V =}1 r 𝜕 𝜕𝑟(𝑟 𝜕V 𝜕𝑟) + 1 𝑟2 𝜕2_V 𝜕𝜗2+ 𝜕2_𝑉 𝜕2_z= 0 若系統為軸對稱，在𝜗分量沒有變化，可以化簡為： ∇2_{V =}𝜕2V 𝜕2_r+ 1 r 𝜕V 𝜕𝑟+ 𝜕2_𝑉 𝜕2_z = 0 --式 2 (二) 電容 電容只和形狀有關，與電壓、電荷量無關（電容定義在此略過不解釋）。 C = 𝑞 𝑉 當電容器的形狀為規則的幾何圖形（如：平板、圓柱、球體），可以使用高斯定律的積 分形式算出 q，以電場內的路徑積分求出 V，便可得到電容值 C。 𝑞𝑒𝑛𝑐 = 𝜀0∯ 𝐸⃗ ∙ 𝑑𝑎 V =ΔU q = − Δ𝑊_𝐸 q = − ∫ 𝑞𝐸⃗ ∙ 𝑑𝑟 但是，推廣到任意形狀時，在數學上就變得不方便計算。於是，可用能量的觀點切入。 當我們把小小的正電荷dq 從電位 Vi移到Vf，我們需要對它做功： 𝑑𝑊 = 𝑈𝑓− 𝑈𝑖 = 𝑉 𝑑𝑞 = 𝑞 𝐶𝑑𝑞 所以，把正電荷q 從電位 Vi移到Vf所需的能量，也就是電容所儲存的位能： 𝑈 = ∫𝑞 𝐶𝑑𝑞 = 𝑞2 2𝐶 = 1 2𝐶𝑉2 我們又知道，電容器中的能量密度是： 𝑢 =𝑑𝑈 𝑑𝜏 = 1 2𝜀0E2 連結兩條式子，可得： 𝐶 =2 ∭ 1 2𝜀0E2𝑑𝜏 𝑉2 = 𝜀0 𝑉2∭ E2𝑑𝜏 --式 3 本程式模擬的目標，便是以式 1 或式 2，加上及式 3，求得電容器之電容值 C。

參、 應用

一、驗證模擬

(一) 平行電板 設計兩平行電容板，兩板均為0.1m×0.1m 之正方形，兩板間距 0.01（m），下盤電位為 5000（V），故此均勻電場 E=5×105_{（V/m）向上。} 1. 理論計算 令高斯面為包住下盤之長方體。 C =𝑄 𝑉 = 𝜀₀∯ 𝐸⃗ ∙ 𝑑𝑎 5 × 103 = 𝜀₀(5 × 105_)0.12 5 × 103 =8.854187818× 10−12（F） 2. 數值模擬 在此可只處理二維的情形，因系統的每個截面形狀均相同。由式1： 𝜕2_V 𝜕2_x+ 𝜕2_V 𝜕2_y = 0 用前面提過的概念，將控制方程式拆為差分式： 𝜕2_V 𝜕2_x+ 𝜕2_V 𝜕2_y = 𝑉_𝑖+1,𝑗−2𝑉_𝑖,𝑗+ 𝑉_{𝑖−1,𝑗} ∆x2 + 𝑉_𝑖,𝑗+1−2𝑉_𝑖,𝑗+ 𝑉_{𝑖,𝑗−1} ∆y2 = 0 加上虛擬時間（注意，我將座標為(i,j)處虛擬時間設為 N，而(i,j)處周圍的點設為 N-1，這就好像拿一個點周圍已知數值的點來算中央這個我們還不知道的點）： 𝑉_𝑖+1,𝑗𝑁−1_{− 2𝑉} 𝑖,𝑗𝑁+ 𝑉𝑖−1,𝑗𝑁−1 ∆x2 + 𝑉_𝑖,𝑗+1𝑁−1_{− 𝑉} 𝑖,𝑗𝑁+ 𝑉𝑖,𝑗−1𝑁−1 ∆y2 = 0 𝑉_𝑖,𝑗𝑁 ₌(𝑉𝑖+1,𝑗𝑁−1+ 𝑉𝑖−1,𝑗𝑁−1)∆y2+ (𝑉𝑖,𝑗+1𝑁−1+ 𝑉𝑖,𝑗−1𝑁−1)∆x2 2(∆x2_{+ ∆y}2₎ 使用混合型邊界條件，上板電位設為1.0（V），下板電位設為5001.0（V），左右兩側 電場水平分量為0，得（備註：將所有的電位值加 1.0，是在不改變電場的狀況下， 方便程式中誤差的討論，見附錄）： 圖2：平行電板電容模擬圖 註：（1）模擬求出之電容值為8.85419 × 10−12_{（F），與理論值之誤差：3.67 × 10}−8_{%。（2）當每處} 電位變化對電位的比值皆小於10−8_{，疊代便停止。}

我們可以把下盤想成是高處、上盤想成平地，方程式所做的事情就只是由我們給定的 「這個區域的輪廓」，將其餘的部分算出來，得到一塊「斜坡」。視覺化後如下： 圖3：由下盤至上盤之電位（左）及電場（右） 註：平行電板間的電場為均勻電場，計算後可知確實如此。 (二) 圓柱形電容器 設計半徑𝑅₁=0.01（m）及𝑅₂=0.02（m）之同軸導體圓柱，長度為 0.05（m），內圓柱電 位為5000（V），外圓柱電位為 0（V）。 1. 理論計算 由高斯定理求出圓柱導體產生的電場。取和此圓柱同軸之圓柱，半徑略大於此圓柱。 V = − ∫ 𝐸⃗ ∙ 𝑑𝑠 𝑅1 𝑅2 = − (_4𝜋𝜀𝑄 0𝐿∫ 1 𝑟𝑑𝑟 𝑅1 𝑅2 ) = − 𝑄 4𝜋𝜀₀𝑙𝑛 |𝑟|| 𝑅1 𝑅₂ C =𝑄 𝑉 = 𝑄 𝑄 4𝜋𝜀₀𝑙𝑛𝑅𝑅2₁ = 4.013037 × 10−12_（F） 2. 數值模擬 此系統可以從兩個方向切入，使之簡化為二維的問題，減少計算量。第一，可以將圓 柱過軸縱切，如此一來，因為圓柱旋轉對稱，只需處理直線的邊界條件，不過方式是 須改為式2 的圓柱座標形式。第二，我們可以將圓柱過軸橫切，對於任一高度而言皆 是對稱的，如此一來，便可向平行電板一樣，用直角坐標的方程，並且要給出圓形的 邊界條件。為了示範圓形的邊界條件，為後續的模擬做準備，在此我用第二種作法。 使用狄利克雷邊界條件，內圓柱電位設為 5001.0（V），外圓柱電位設為 1.0（V）， 得：

圖4：圓柱形電容器之電容模擬圖 註：（1）模擬求出之電容值為4.00936 × 10−12_{（F），與理論值之誤差：9.16 × 10}−2_{%。我們可以發} 現到，圓柱形電容器的數值模擬誤差較平行電板的數值模擬誤差多了6 個數量級。此現象很可能是因 為網格的設置為方格，當邊界條件為圓形的時，誤差相對來說較大。若要求更精準的計算，可採用不 同形狀的網格（在此不深入解釋）。（2）當每處電位變化對電位的比值皆小於10−8_{，疊代便停止。} 圖5：由內圓柱至外圓柱之電位（左）及電場值（右） 註：就像平行電板，由高電位到低電位像是一個斜坡。任意取過圓心、內到外的路徑來觀察，可發現 在內圓柱內，因無電位差，電位像是一個台地，電場為零。而在兩圓柱間的電位變化(電場值)，由內而 外是由陡峭逐漸平緩，猶如滑雪斜坡坡道一般。

(三) 球形電容器 設計半徑𝑅₁=0.016（m）及𝑅₂=0.02（m）之同心導體球殼，小球殼電位為 5000（V）， 大球殼電位為0（V）。 1. 理論計算 由球殼定理，小球殼產生的電場可簡化為球心之點電荷發出的電場。計算內外球殼電 位差時，選定積分路徑為垂直於內球殼，由外球殼到內球殼： V = − ∫ 𝐸⃗ ∙ 𝑑𝑠 𝑅1 𝑅2 = − (_4𝜋𝜀𝑄 0∫ 1 𝑟2𝑑𝑟 𝑅1 𝑅2 ) = − 𝑄 4𝜋𝜀₀(− 1 𝑟)| 𝑅1 𝑅₂ C =𝑄 𝑉 = 𝑄 𝑄 4𝜋𝜀₀(𝑅1₁−𝑅1₂) = 8.90120 × 10−12_（F） 2. 數值模擬 雖然此系統看似是三維的問題，但其實邊界條件對於任何一個通過球心的軸都是對 稱的。我們可以任取過球心的直線為軸，以圓柱座標來解決這個問題。為方便，取過 球心的鉛直線為軸。邊界條件完全和圓柱形電容相同，只是方程式換了一個樣子。由 式2： 𝜕2_V 𝜕2_r+ 1 r 𝜕V 𝜕𝑟+ 𝜕2_𝑉 𝜕2_z= 0 改 r 為 x，設𝑥0為所取之軸的 x 座標，以 i 表示 x(水平軸)，j 表示 z(垂直軸)。 方程式的第二項原本應近似為： 1 r 𝜕V 𝜕𝑟 ≈ { 1 |x − 𝑥0| 𝑉_{𝑖,𝑗−1}− 𝑉_𝑖,𝑗+1 2∆x , x < 𝑥0 1 |x − 𝑥0| 𝑉_𝑖,𝑗+1− 𝑉_{𝑖,𝑗−1} 2∆x , x > 𝑥0 為求簡潔與方便（負負得正），我將它寫成如下形式（注意，程式中必須刻意將格子 點與𝑥₀座標錯開，否則|x − 𝑥₀| = 0，會出現 0 在分母的情形）： 𝑉_𝑖+1,𝑗−2𝑉_𝑖,𝑗 + 𝑉_{𝑖−1,𝑗} ∆x2 + 1 x − 𝑥₀ 𝑉_𝑖,𝑗+1− 𝑉_{𝑖,𝑗−1} 2∆x + 𝑉_𝑖,𝑗+1−2𝑉_𝑖,𝑗+ 𝑉_{𝑖,𝑗−1} ∆z2 = 0 移項，並加上虛擬時間，得： 𝑉_𝑖,𝑗𝑁 ₌(𝑉𝑖+1,𝑗 𝑁−1_{+ 𝑉} 𝑖−1,𝑗𝑁−1)∆z2+ (𝑉𝑖,𝑗+1𝑁−1+ 𝑉𝑖,𝑗−1𝑁−1)∆x2+∆x 2_∆z2 x − 𝑥₀ (𝑉_𝑖+1,𝑗𝑁−1 + 𝑉_{𝑖−1,𝑗}𝑁−1₎ 2∆z 2(∆x2_{+ ∆z}2₎

其中，𝑥₀為所取之軸的 x 座標。 𝑉_𝑖,𝑗𝑁 ₌(𝑉𝑖+1,𝑗 𝑁−1_{+ 𝑉} 𝑖−1,𝑗𝑁−1)∆z2+ (𝑉𝑖,𝑗+1𝑁−1+ 𝑉𝑖,𝑗−1𝑁−1)∆x2+∆x 2_∆z2 x − 𝑥₀ (𝑉_𝑖+1,𝑗𝑁−1 + 𝑉_{𝑖−1,𝑗}𝑁−1₎ 2∆x 2(∆x2_{+ ∆z}2₎ 使用狄利克雷邊界條件，內圓柱電位設為 5001.0（V），外圓柱電位設為 1.0（V）， 得： 圖6：球形電容器之電容模擬圖 註：模擬求出之電容值為8.90925 × 10−12_{（F），與理論值之誤差：9.04 × 10}−2_{%。我們注意到，球形} 電容器的數值模擬誤差也較平行電板的數值模擬誤差多了 6 個數量級，且甚至比圓柱形電容器的誤差 還要大一些。這是能夠理解得，因為圓柱形電容器只有在切分圓形的部分不準，在垂直方向（即，圓 柱座標的 Z 軸）是準確的；而球型電容器在三個維度都是不準的，誤差然就更大了。同樣地，若設置 更好的網格，此問題能獲改善。

圖7：球形電容器之電容模擬之收斂情形 註：此圖為疊代的過程中，每十次（即附錄程式碼中的 gap）計算一次電容，觀察其收斂情形。藍色虛 線為模擬的電容值，而紅色實線為理論值。由於誤差較前兩例大，因此將收斂條件變得更嚴格：當每 處電位變化對電位的比值皆小於10−11_{，疊代便停止。} 圖8：由內球至外球之電位（左）及電場（右） 註：就像先前舉的兩個例子，由高電位至低電位（任意取過球心、內到外的路徑來觀察），也像是一 個斜坡。

二、求解不規則形電容

由平行電容板公式（C =𝜀0A d ，A 為平行板面積、d 為兩者間距）可知，當兩盤越靠近， 電容就越大。那麼，若是只有部分突起變形呢？電容值是否也會上升？此難以由理論 計算的狀況，可以先嘗試粗略地估算；若必須求出精確的值，可由程式得到答案。 設計以下形狀：平板上有一個半橢球形金屬球殼（等電位體）突起，兩盤均為直徑0.1 （m）之圓形，橢球中心位於兩盤上，兩圓盤間距 0.01（m），下盤電位為 0（V），上盤 電位為5000（V）。令半橢球形在下盤上的軸不變，為 0.01（m），而它的高（垂直於下 盤的半軸長）由0.002 變化至 0.009（m），討論逐漸變化時，電容值的變化。

(一) 定性討論 在定性討論時，可以嘗試將半橢球凸起部分假想為另一平板，如下圖： 圖9：定性想法概念圖 註：其中 d1> d2，而其中 A2為半橢球在下盤所占的面積，A1為下盤剩餘的面積；(b)相當於三個小電容 的並聯，其電容值為三個電容值相加。 經果簡單的計算，我們可以知道，(a)之電容值<(b)之電容值<(c)之電容值，因為： 𝜀₀(𝐴₁+ 𝐴₂) 𝑑₁ < 𝜀₀𝐴₁ 𝑑₁ + 𝜀₀𝐴₂ 𝑑₂ < 𝜀₀(𝐴₁+ 𝐴₂) 𝑑₂ 故大略得知，此種形狀之電容應介於「間距0.01(m)之平行電板電容」與「間距和半橢 球頂端與上盤距相同之平行電板電容」，且隨著半橢球離上盤越來越近，電容值會越來 越大。 (二) 程式模擬 每0.001（m）計算一次。使用上述之圓柱座標，使用混合型邊界條件，下盤設為 5001.0 （V），上盤設為 1.0（V）。此僅畫出當中的兩個呈現： 圖10：電容模擬圖 註：上圖，半橢球之半短軸（垂直）為 0.002（m）；下圖，半橢球之半長軸（垂直）為 0.009（m）。

圖11：電容變化圖 註：定義橫坐標 h 代表半橢球上方頂點與上盤的距離。隨著距離縮短，電容值逐漸上升。 為了進一步知道電容與凸起形狀的關係，思考因平行電板電容正比於兩盤間距的倒數 ，又平板只有一小塊形狀改變，大部分區域的電場不受影響，也許會想延伸定性討論中 的想法（見圖9(b)），將凸起部分視為一被提高的平板，用C(h) − 𝐶0=𝐾_h來擬合這些數 據點；但是，考慮到若只用距離 h 不能代表凸起部分與上盤的距離，不妨求出凸起部 分與上盤之平均距離ℎ̅，再嘗試以C(ℎ̅) − 𝐶0 =𝐾_ℎ̅來擬合此數據。 平均距離ℎ̅可由上盤與半橢球之間的體積除以半橢球在下盤所占的面積來求得，示意圖 如下圖： 圖12：平均高度計算示意圖

圖13：電容值擬合圖 註：藍色點為程式模擬數據點，紅線是由擬合軟體 Logger Pro 計算得出，待擬合參數為 C0及 K，橫坐 標 h-bar (ℎ̅)代表半橢球上方頂點與上盤的平均距離。 擬合結果如下： 𝐶(ℎ̅) = 6.682 × 10−12₊2.404 × 10−15 ℎ̅ 同時，因擬合想法是由定性討論的模型出發，我們可以簡單地驗證擬合結果：若只考慮 扣除凸起部分的平行電板（如圖9(b)中的 A1），其電容值為ε0_dA1 1 = 6.885 × 10 −12_，與擬 合得出的6.682×10-12_{相近，可知程式模擬結果及擬合結果是合理的。可以將}2.404×10−15 ℎ̅ 理解成ε0_𝛼ℎ̅A2（其中α 為一參數，表示ℎ̅與等效平行板電容之間距，αℎ̅對應圖 9(b)中的 d2 ）。解得α=0.2893，意味著能將中央凸起比擬為間距為 0.2893ℎ̅的平行電板。 (a) (b) 圖14：定量分析概念圖 註：(a)分析第一步驟，將 h 變換為平均距離ℎ̅，進行擬合；(b)將電容凸起部分類比為距離相距 0.2893ℎ̅ 的平行電板電容。(a)和(b)兩形狀之電容值相等。

0.2983 這個值看似有點小，但其實是可以理解的。因為半橢球頂端的曲率半徑相較平 板來說很小，單位面積內能累積更多的電荷，電場也就更強（類似尖端放電的概念）， 由式3 可知 C∝∭ E2_{𝑑𝜏，便可以大略猜到半橢球凸起讓電容值上升很快的現象。.} 經由定性和定量的分析，不但知道電容值隨著突起的半橢球越來越大而漸增，得出了 電容的精確值，還得到擬合結果，可以簡易地預測其他半橢球凸起形狀的電容值，解決 了原先的疑惑。

肆、 結語

本程式的緣起，是由於姜理元同學在專題研究上有求解電容的需要。我剛好有一些解拉 普拉斯方程的經驗，經過一些嘗試，成功的解決了他在理論模型數學計算上的困難。透過模 擬，不但可將物理現象視覺化，讓人易於了解，甚至可以經由數值分析求得未知的物理量。 物理當中，程式模擬幾乎是必備的工具，希望這篇文章對讀者有所幫助。

致謝

感謝中研院博士後研究員鄭明宏博士在程式語言上的指導、姜理元同學在電學方面提供 的諮詢及賴奕帆老師的建議。

參考文獻

1. 徐力宏、方淳民、黃仁偉（譯）（民100）物理（下）（第九版）（原作者：David, H., Robert, R, Jearl, W.）。新北市：全華圖書。

2. Erwin Kreyszig. (1983). Advanced Engineering Mathematics(5e). Taipei: The Tan Chiang Book Company.

3. Griffiths, David J.(1999). Introduction to Electrodynamics. New Jersey: Prentice-Hall, Inc. 4. Lorena A. Barbarra. (2013). 12steps to Navier-Stokes. retrieved from：

http://nbviewer.jupyter.org/github/barbagroup/CFDPython/tree/master/lessons/ 5. Matplotlib. (2017). Retrieved from：https://matplotlib.org/

附錄

一、平行電板電容模擬程式碼

'''1. 輸入函式庫''' #--- import sys reload(sys) sys.setdefaultencoding('utf8') import numpy as np import matplotlib.pyplot as plt from sympy import init_printing init_printing() import time tStart = time.time() #--- '''2. 定義參數和矩陣''' #--- epsilon0 = 8.854187817620e-12 L = 0.01 #平板高度 W = 0.1 #平板寬度 V0 = 5000.0 #高電位 nx, ny = 800, 80 #網格 x,y 方向之格線數 dx = W/(nx-1) #網格邊長(寬) dy = L/(ny-1) #網格邊長(直) x = np.linspace(0, W, nx) #網格 x 座標 y = np.linspace(0, L, ny) #網格 y 座標 X, Y = np.meshgrid(x, y) #生成座標矩陣(畫電位用) phi = np.full((ny,nx), 1.0) #儲存電位的陣列 ###生成座標矩陣(畫電場用)，其座標都在電位座標格子點的中央 X1 = (X[0:ny-1,0:nx-1] + X[0:ny-1,1:nx] + X[1:ny,0:nx-1] + X[1:ny,1:nx])/4 Y1 = (Y[0:ny-1,0:nx-1] + Y[0:ny-1,1:nx] + Y[1:ny,0:nx-1] + Y[1:ny,1:nx])/4 #--- '''3. 邊界條件''' #--- #將所有的電位值加 1.0，在不改變結果的狀況下，方便誤差的討論(避免 0 在分母) phi[0, :] = V0+1.0 #上邊界(在圖中上下顛倒) phi[ny-1, :] = 0.0+1.0 #下邊界(在圖中上下顛倒) phi[1:ny-1, 1] = phi[1:ny-1, 0] + (0.0)*(1.0/100) #左邊界 phi[1:ny-1, nx-2] = phi[1:ny-1, nx-1] - (0.0)*(1.0/100) #右邊界 #--- '''4. 主程式碼''' #--- count = 0 #計數疊代多少次 C = [] #儲存每次疊代的電容值 ratio_container = [] #儲存「一次疊代前後的差與疊代前的電位比值」

times = int(1e2) #疊代次數上限 for N in range(1, times):

pre_phi = np.copy(phi) #複製上一次疊代的值 #拉普拉斯方程(直角坐標)

phi[1:ny-1, 1:nx-1] = \

(dx**2*(pre_phi[2:ny, 1:nx-1]+pre_phi[0:ny-2, 1:nx-1])+\

dy**2*(pre_phi[1:ny-1, 2:nx]+pre_phi[1:ny-1, 0:nx-2]))/(2*(dx**2+dy**2)) count += 1 '''邊界條件''' phi[0, :] = V0+1.0 #上邊界(在圖中上下顛倒) phi[ny-1, :] = 0.0+1.0 #下邊界(在圖中上下顛倒) phi[1:ny-1, 1] = phi[1:ny-1, 0] + (0.0)*(1.0/100) #左邊界 phi[1:ny-1, nx-2] = phi[1:ny-1, nx-1] - (0.0)*(1.0/100) #右邊界 '''電場 = -gradient(電位)，為了方便，另 U=Ex而 V=Ey'''

U = -1*((phi[0:ny-1,1:nx]+phi[1:ny,1:nx])-(phi[0:ny-1, 0:nx-1]+phi[1:ny, 0:nx-1]))/(2*dx) V = -1*((phi[1:ny, 0:nx-1]+phi[1:ny, 1:nx]) - (phi[0:ny-1, 0:nx-1]+phi[0:ny-1, 1:nx]))/(2*dy) inte = []

'''E**2 對 2*pi*r 和 z 積分'''

integration = (U**2 + V**2)*(dx*dy*W) C.append(epsilon0/V0**2*np.sum(integration))

'''檢驗電位是否收斂(每處的電位值在疊代前後的差與迭代前的電位比值小於 err)''' dphi = phi - pre_phi

ratio = dphi[1:ny-1,0:nx]/phi[1:ny-1,0:nx]

ratio_container.append(np.max(dphi[1:ny-1,1:nx-1]/phi[1:ny-1,1:nx-1])) err = 1e-8#當比值小於 err，停止運算

if np.all(ratio <= 1e-8): print 'converges' break print count #--- '''5. 繪圖''' #--- ###分層設色圖

fig = plt.figure(figsize=(20,5), dpi=200) #建立「繪圖板」 plt.axis('Scaled') #使 x,y 座標單位相同

extent = [0, W, 0, L] #設定繪圖範圍

magnitude = 10**np.floor(np.log10(V0)) #電位大概的數量級

levels = np.arange(0.0, V0/magnitude , (V0/magnitude)/20) #分層設色圖的間隔 cmap = plt.cm.get_cmap("jet") #色系

cmap.set_under(color=(0,0,0.4)) #若值超過分層設色圖的上限，將顏色訂為(R,G,B) cont = plt.contourf(X, Y, phi/magnitude\

,levels=levels, cmap=cmap, linewidths=1, extend='both')

#CS = plt.contour(X, Y, phi/magnitude, linewidths=1) #Equal planes of electric potential #等位線， 可畫可不畫

cbar = plt.colorbar(cont)

cbar.set_label('Electric potential $\mathregular{(10^{%r} Volt)}$'%int(np.log10(magnitude)), fontsize=18, rotation=270, labelpad=30) #圖表基本資料

plt.imshow(phi, extent=extent, origin="lower") #畫圖，lower 表示「冷色系為值小者，暖色系為 值大者」

###電場 = -gradient(電位)

U = -1*((phi[0:ny-1, 1:nx]+phi[1:ny, 1:nx]) - (phi[0:ny-1, 0:nx-1]+phi[1:ny, 0:nx-1]))/(2*dx) V = -1*((phi[1:ny, 0:nx-1]+phi[1:ny, 1:nx]) - (phi[0:ny-1, 0:nx-1]+phi[0:ny-1, 1:nx]))/(2*dy) magnitude_UV = 2*10**np.ceil(np.log10(np.sum(V)/(nx*ny))) #電場大概的數量級

interval = 9 #每 interval 個資料點畫一個箭頭表示電場

Q=plt.quiver(X1[0:ny:interval, 0:nx:interval], Y1[0:ny:interval, 0:nx:interval],\

U[0:ny:interval, 0:nx:interval], V[0:ny:interval, 0:nx:interval], scale = magnitude_UV, scale_units='inches', width = 0.0015)

#畫箭頭圖例

qk = plt.quiverkey(Q, 0.9, 1.08, magnitude_UV, r'$\mathregular{%r\times10^{%r}V/m}$' %(magnitude_UV/10**int(np.log10(magnitude_UV)),int(np.log10(magnitude_UV))),\ labelpos='E', fontproperties={'weight': 'bold', 'size':'large'})

#畫座標

plt.axis([0, W, 0, L]) #文字樣式

label_font = {'fontsize':16, 'verticalalignment':'top','horizontalalignment' :'center', 'verticalalignment' :'center'}

#x,y 軸、圖表標題

plt.xlabel('$(m)$', label_font) plt.ylabel('$(m)$', label_font) plt.title('electric field',fontsize=20)

fig.savefig('parralel dsik(error=1e-8)%rtimes%r, %rV, value=%r.png'%(W,L,V0, C)) #存檔 #--- '''6. 誤差討論''' #--- ratio_container = np.asarray(ratio_container) fig_ratio = plt.figure() line, = plt.semilogy(np.linspace(0,count,count),ratio_container,'b-',linewidth=2.0) plt.xlabel('times',label_font) plt.ylabel('ratio',label_font) plt.show() fig_ratio.savefig('ratio.png') #存檔 #---

plt.show()

tEnd = time.time()

print "total time is", tEnd - tStart #算出程式運算時間

二、

球形電容器模擬程式碼（僅列出與平行電板電容模擬程式碼不同處） (一)'''2. 定義變數和矩陣''' 1. 宣告新的參數 R = 0.02 #外球半徑 r = 0.016 #外球半徑 2. nx 刻意不用漂亮的數字，使座標轉換時，圓柱座標的垂直軸不與原本的網格線重疊， 避免在疊代中出現 0 在分母的情形 nx, ny = 404, 401 #the size of the grid

3. 加上用於座標轉換的陣列

###座標平移用之陣列（轉換為圓柱座標） calibration_origin = np.full((ny, nx), center[0])

(二)'''3. 邊界條件'''改為'''3-1 設定邊界位置'''和'''3-2 邊界條件''' '''3-1 設定邊界位置'''

###定義計算範圍(x 座標)

extent = [range(int((center[0]-r)/dx), int(np.ceil(((center[0]+r)/dx))+1)),

range(int((center_outer[0]-R)/dx), int(np.ceil((center_outer[0]+R)/dx)+1))] #left and right edge of the inner electrode #left and right edge of the outer electrode

extent = np.asarray(extent)

###(pos_x, pos_y)為邊界條件的位置 pos_x = np.linspace(0.0 ,W, nx)

pos_y = np.full((2,nx,2),center[1]) #pos_y[0,:,:]表示外部邊界，pos_y[1,:,:]表示外部邊界 def positive(x):#定義一函數限制 pos_y 的範圍

if x >= 0: return x if x < 0: return 0.0 #由圓的一般式:(x-center[0])**2+(y-center[1])**2=r**2 計算 pos_y for i in extent[0]:

pos_y[0][i][1] = center[1] + np.sqrt(positive(r**2 - (pos_x[i]-center[0])**2)) pos_y[0][i][0] = center[1] - np.sqrt(positive(r**2 - (pos_x[i]-center[0])**2)) for i in extent[1]:

pos_y[1][i][1] = center_outer[1] + np.sqrt(positive(R**2 - (pos_x[i]-center_outer[0])**2)) pos_y[1][i][0] = center_outer[1] - np.sqrt(positive(R**2 - (pos_x[i]-center_outer[0])**2)) #定義一陣列表示電極間空間的 y 座標區間

vertical_extent_grid = [map(lambda x: int(round(x)), pos_y[0][:, 1]/dy), map(lambda x: int(round(x)), pos_y[1][:, 1]/dy)]

#(依照上一部分 pos_y 的求法，pos_x→pos_y 為一對一，使旁邊出現空隙出現)

def pos_boundary(l ,m): #參數'l'=0 或 1，分別表示內球或外球的 x 座標範圍；參數'm'=0 或 1，分別表示上半球或下半球

pos_bound = [] ext = extent[l]

vertical_extent = map(lambda x: int(round(x)), pos_y[l][ext, m]/dy) pos_bound.append([ext[0], vertical_extent[0]])

for i in range(0,len(ext)/2):#左半部

if vertical_extent[i+1] == vertical_extent[i]: pos_bound.append([ext[i+1], vertical_extent[i+1]]) else:

sign = (vertical_extent[i] - vertical_extent[i+1])/abs(vertical_extent[i] - vertical_extent[i+1]) gap_grid = abs(vertical_extent[i] - vertical_extent[i+1])

for j in range(gap_grid):

pos_bound.append([ext[i+1], vertical_extent[i]-sign*(j+1)]) j+=1

for i in range(len(ext)/2, len(ext)-1): #右半部 if vertical_extent[i] == vertical_extent[i+1]: pos_bound.append([ext[i], vertical_extent[i]]) else:

sign = (vertical_extent[i] - vertical_extent[i+1])/abs(vertical_extent[i] - vertical_extent[i+1]) gap_grid = abs(vertical_extent[i] - vertical_extent[i+1])

for j in range(gap_grid): pos_bound.append([ext[i], vertical_extent[i]-sign*j]) j+=1 pos_bound.append([ext[len(ext)-1], vertical_extent[len(ext)-1]]) return np.asarray(pos_bound) #四個二維陣列表示個半球的邊界條件位置，第一行為 x 座標，第二行為 y 座標 pos_boundary00 = pos_boundary(0,0) pos_boundary01 = pos_boundary(0,1) pos_boundary10 = pos_boundary(1,0) pos_boundary11 = pos_boundary(1,1) '''3-2 邊界條件''' #將所有的電位值加 1.0，在不改變結果的狀況下，方便誤差的討論(避免 0 在分母)，見'''4. 主程式碼'''的'''檢驗電位是否收斂''' phi[pos_boundary00[:,1], pos_boundary00[:,0]] = V0+1 phi[pos_boundary01[:,1], pos_boundary01[:,0]] = V0+1 phi[pos_boundary10[:,1], pos_boundary10[:,0]] = 0.0+1 phi[pos_boundary11[:,1], pos_boundary11[:,0]] = 0.0+1 phi[0, :] = 0.0+1 #上邊界(在圖中上下顛倒) phi[ny-1, :] = 0.0+1 #下邊界(在圖中上下顛倒) phi[:, 0] = 0.0+1 #左邊界 phi[:, nx-1] = 0.0 +1 #右邊界

(三)'''4. 主程式碼'''的球形電容器版本 '''4. 主程式碼'''

count = 0

C = [] #儲存每次疊代的電容值

ratio_container = [] #儲存「疊代前後 phi 的差與疊代前的 phi 比值」的最大值 gap = 10 #每 gap 次檢驗一次收斂情形

times = int(1e2) #疊代次數上限 for N in range(1, times):

pre_phi = np.copy(phi) #an array contains phi in the previous iteration ###拉普拉斯方程(圓柱坐標)

phi[1:ny-1, 1:nx-1] = \

((pre_phi[1:ny-1, 2:nx]+pre_phi[1:ny-1, 0:nx-2])*(dy**2) \ +(pre_phi[2:ny, 1:nx-1]+pre_phi[0:ny-2, 1:nx-1])*(dx**2) \ - dx*dy**2/(X[1:ny-1, 1:nx-1] \ - calibration_origin[1:ny-1, 1:nx-1])/2*(pre_phi[1:ny-1, 2:nx]-\ pre_phi[1:ny-1, 0:nx-2])/2)/(2.0*(dx**2+dy**2)) count += 1 '''邊界條件''' phi[pos_boundary00[:,1], pos_boundary00[:,0]] = V0+1 phi[pos_boundary01[:,1], pos_boundary01[:,0]] = V0+1 phi[pos_boundary10[:,1], pos_boundary10[:,0]] = 0.0+1 phi[pos_boundary11[:,1], pos_boundary11[:,0]] = 0.0+1 phi[0, :] = 0.0+1 #up phi[ny-1, :] = 0.0+1 #down(enclosed) phi[:, 0] = 0.0+1 phi[:, nx-1] = 0.0+1

#if I didn't adjust all the boundaries by 1.0, division by 0.0 occurs while computing error. if N%gap == 0:

span_container = [] '''電場 = -gradient(電位)'''

U = -1*((phi[0:ny-1, 1:nx]+phi[1:ny, 1:nx]) - (phi[0:ny-1, 0:nx-1]+phi[1:ny, 0:nx-1]))/(2*dx) V = -1*((phi[1:ny, 0:nx-1]+phi[1:ny, 1:nx]) - (phi[0:ny-1, 0:nx-1]+phi[0:ny-1, 1:nx]))/(2*dy) '''E**2 對 2*pi*r 和 z 積分''' inte = [] for i in extent[1][0:len(extent[1][:])-1]: #外球的 x 座標範圍-1，表示電場的 x 座標索引值 if vertical_extent_grid[0][i] == vertical_extent_grid[1][i]: continue inte.append( sum((U[vertical_extent_grid[0][i]:vertical_extent_grid[1][i], i]**2 + V[vertical_extent_grid[0][i]:vertical_extent_grid[1][i], i]**2)*(dx*dy)*\

(X1[vertical_extent_grid[0][i]:vertical_extent_grid[1][i], i]- \ np.full((vertical_extent_grid[1][i]-vertical_extent_grid[0][i]), center[0])))) ###instant look of C C.append(2*epsilon0/V0**2*(np.pi*sum(np.absolute(inte)))) '''檢驗電位是否收斂(每處的電位值在疊代前後的差與迭代前的電位比值小於 err)''' for i in extent[1]: #外球的 x 座標範圍 if vertical_extent_grid[0][i] == vertical_extent_grid[1][i]: continue

dphi = phi - pre_phi

span_container.append(np.asarray(dphi[vertical_extent_grid[0][i]:vertical_extent_grid[1][i], i])/ \ np.asarray(phi[vertical_extent_grid[0][i]:vertical_extent_grid[1][i], i])) ratio1 = [] #儲存每一直行中 dphi/phi 的最大值 for i in range(len(span_container)): ratio1.append(max(span_container[i-1])) ratio_container.append(max(ratio1)) err = 1e-11 #當比值小於 err，停止運算 if ratio_container[count/gap-1] <= err: print 'converges'

tEnd = time.time()

print "total time is", tEnd - tStart print count

Solving the Capacitance with various boundary conditions

through python

Chih-Chieh Wu

*

, Li-Yuan Chiang, and Yi-Fan Lai

Taipei Municipal Jianguo Senior High School *Corresponding author: jackwu9923@gmail.com

Abstract

This essay briefly introduces finite difference method, which was adopted to solve a Laplace’s equation of electric potential with a view to promoting the simple but useful simulating technique. The basic concept of finite difference method is first explained, namely, approximation through Taylor expansion, the setting of mesh, and the three form of boundary conditions. Next, comparisons are made between theoretical and simulated capacitance of parallel-plate capacitor, cylindrical capacitor, and spherical capacitor. Last, the capacitances of non-standard-shape capacitors are computed, demonstrating the usefulness of the method.