有限元素法與Poisson Boltzmann方程式驗證

3-1 有限元素法簡介

有限元素法(Finite Element Method, FEM)是使用一不連續之函數來近似任一連續 量。換句話說，經由我們所建立的函數，可以描述一整體結構的行為。有限元素法的操 作方式是將一連續體分割成數個四面體(2D/三角形)、六面體(2D/矩形)或五面體，稱之 為元素（element），而元素之邊界點稱之為節點（node），例如四面體就有四個節點，

六面體則有八個節點，每個節點上皆可用數學方程式來描述，稱之為內插函數方程式

（interpolation equation），藉由這些內插函數方程式表達該連續體之分析行為，此群有 限個方程式之解稱為內插近似解（interposition approximation）。只要連續體之場變數（應 變、應力、溫度、電場、磁場、濃度…等）與相關正確，則近似解可視為精確解。

有限元素法的優點包括：

1. 元素的材料性質不需相同，例如非均質、非等向性，只要設定元素的條件，就 可以完成近似解求解。

2. 不規則形狀的邊界，能用直邊的小元素作近似估計，或用二階函數元素趨近。

3. 元素的大小、疏密，可經由實際所要計算的場來做調整。例如變化劇烈的區域 可以使用較密的元素來填佈。

4. 容易處理混合型的邊界條件。

因此，有限元素法被廣泛地應用在各種領域上，包括了固力、流力、熱力、製造以 及結構等。本研究所採用的是三維的元素網格，更可以實際模擬出立體的場變化。

3-2 電雙層理論與有限元素法的離散

當固、液、氣三相中之任二相彼此接觸時，於界面上因材料之電負度差異，會產生 分子間的電荷分佈不均，將使界面材料表面呈現帶電荷或極化情形。此材料界面上極性 會吸引帶異電性的電荷吸附與聚集，以形成所謂的電雙層(Electric double layer, EDL)。

根據其理論：電勢(electrical potential) ψ 與 每單位體積內的靜電荷密度ρ_e 的三維普 松方程式(Poisson’s equation)可以描述為

0 2

2

2 2

2 2 2

εε ψ ρ

ψ

ψ ψ ^e

z y

x

=−

∂ +∂

∂ +∂

∂

=∂

∇ (3-1)

其中ε 是介電常數(dielectric constant of the solution)，ε₀ 則是真空的介電常數 (permittivity of vacuum)。而離子的濃度呈現對稱的形式：

)

將式(3-3)帶入式(3-1)，則可以得到非線性的 Poisson-Boltzmann 方程式 )

接著，我們定義Debye-Huckel parameter

T

加入四面體元素(tetrahedral element, 如圖 3-1 所示)，一階 Shape function 可以被展 開成為

)

將式(3-10)至式(3-14)帶入式(3-9)可得

應用Gauss divergence 定律，上式成為

∑∫∫∫

最後，我們使用式(3-7) trial solution 帶入上式，可得

∑∫∫∫

若將

∫∫∫ ^dxdydz

^=V(e)帶入式(3-20)，則可得

最後我們得到Poisson-Boltzmann 方程式之有限元素離散的矩陣形式為

{ } { }

^{( ) ( )} 成使用二階Shape function，則可以使用較大的元素，或者以原來的元素量，得到更精準 的結果。圖3-2 則是二階有限元素，可以跟圖 3-1 的一階有限元素做比較。

⎪ ⎪

成一階函數與二階函數共用shape function。例如：

L

j

N

_i

L

_j

N

_i

4 ) 2 1 4 9 (

) (

9

) 2 2 )( 1 (

1 2 2 1

2

* 1 2

* 1 2

* 1

2

2 1 2

* 1 2

* 1 2

* 1 1 1

+ + −

= +

+

− +

= +

∇

⋅

∇

L v L

d c b

v

L d

c φ b

φ

經由積分公式

d V c b a

d c b dv a

L L L

V

L

d c b

a 6

)!

3 (

!

!

!

!

4 3 2

1 = + + + +

∫

就可以得到

∫∫∫

^{∇φ1⋅∇φ1}

^dxdydz

⁼⁽

^b

^1*2+

^c

₉^1*

_v

^{2 +}

^d

^1*2) ₂₀³

其他項的推導請參閱附錄二的推導。

3-4 Poisson-Boltzmann 方程式驗證

下圖為我們的測試case：球體(radius a=0.325μm)靠近一平板(distance h = 2.5μm)的網格 分佈圖。

圖3-3 表面網格分布圖

接著，我們與先前的論文做比較，證明了我們所撰寫的程式無誤。

圖3-4 Poisson-Boltzmann 方程式驗證 下圖是測試case 的表面電荷分布：

圖3-5 表面電荷分布圖

下表是測試cases，包括了元素總 nodes 數量與四面體元素數量。

Mesh Number of nodes Number of elements 00 4,158 19,958 01 10,846 56,399 02 40,577 219,115 03 50,142 270,766

我們分別以一階、二階shape function 執行 mesh00 至 mesh03 的 case，繪出了圖 3-6。

可以發現若僅以mesh00 的 case 執行二階 shape function，就可以趨近於以 mesh04 執行 一階shape function 的 case。換句話說，只要使用 nodes 數量 4,158(elements 19,958)的二 階函數，就可以得到像使用nodes 數量 50,142(elements 270,766)的一階函數效果，足足 差了十倍左右的nodes 數量與 element 元素數量。

圖3-6 一階與二階 Shape function 比較

我們將實際的平行化執行時間來做比較：nodes 數量 4,158(elements 19,958)的二階 函數，與nodes 數量 50,142(elements 270,766)的一階函數，兩者計算相差不到四秒(二階 函數花103.5 杪、一階函數 99.8 秒)。

下圖是六個CPU 切割網格計算：元素與網格點分布的情形。在 CPU 的邊界必須有 很好的網路溝通與交換效能，各計算主機才能很快地彼此交換計算資料。

圖3-7 Poisson-Boltzmann 方程式的 subdomains

在文檔中 有限元素平行化三維Poisson-Boltzmann Equation程式之應用 (頁 27-37)