誘電率が異なる場合の静電場の数値計算

山本昌志^{∗}

2007 年 7 月 10 日

概要

一部に誘電体を挿入した場合のポアソン方程式の数値計算の方法を示す．

1 方程式

誘電体があっても静電場では，

∇ · D = ρ  ∇ × E = 0  (1)

を満足する．もちろん， D = εE が成り立つ．とくに，図 1 のように一部の領域の誘電率が不連続に変化— 誘電体が連続的に変化しない—場合，

∇ · E = \frac{ρ}{ε_{1}}  ∇ × E = 0  領域 Ω_{1}  (2)
∇ · E = \frac{ρ}{ε_{2}}  ∇ × E = 0  領域 Ω_{2}  (3)

となる．この場合でも， ∇ × E = 0 が成り立つので，電場はスカラーポテンシャル ϕ の勾配と書ける．

E = −∇ϕ  (4)

もっと言えば，誘電体が連続に変化しても，静電場の問題ではポテンシャルの計算に帰着できる．なぜならば，静電場の問題はいつでも ∇ × E = 0 が成り立つからである．

式 (2) と式 (3) を満足させるためには，

∇^{2}ϕ = −\frac{ρ}{ε_{1}}  領域 Ω_{1}  (5)
∇^{2}ϕ = −\frac{ρ}{ε_{2}}  領域 Ω_{2}  (6)

となる必要がある．もし，それぞれの領域に真電荷が無いとすると，ポテンシャルの方程式は，

∇^{2}ϕ = 0  領域 Ω_{1}，領域 Ω_{2}  (7)

∗国立秋田工業高等専門学校 電気情報工学科

となる．ただし，この方程式を解く場合それぞれの領域で計算する必要がある．そして，誘電率が異なる二つの境界でポテンシャルの値を同一にしなくてはならない．これは，それぞれの領域でラプラス方程式を解く問題で，境界条件の設定が重要になる．たぶん，誘電体の境界でポテンシャルを設定することは困難であろう．

そこで，別々にラプラス方程式を解くのではなく，式 (1) とスカラーポテンシャルの式 (4) から，直接得られる

∇ · (ε∇ϕ) = −ρ  (8)

を計算する方がよい．

図 1: ここで取り扱う問題; 誘電体が不連続に変化する領域の数値計算．

2 数値計算方法

次のような条件の静電場を計算する方法について検討を行う．

• 計算領域に真電荷は無い．

ρ = 0  (9)

• 誘電率は不連続に変化する．

このような条件のもと，差分法と有限要素法 (変分法)，有限積分法の数値計算について簡単にコメントする．

2.1 差分法

ここで計算する静電場の方程式は，式 (1) とスカラーポテンシャルから，

∇ · (ε∇ϕ) = 0  (10)

導くことができる．これは 2 階の微分方程式になっており，差分の式はとなりの要素を含めた計算になる．隣の要素まで含めるとなると，一つの式の中に誘電率が異なる部分が生じる (図 2)．この取り扱いは面倒である．

あるいは，それぞれの誘電体中でラプラス方程式

∇^{2}ϕ = 0  (11)

を計算してもスカラーポテンシャルを得ることもできる．しかし，実際この式を計算するこは，簡単ではない．二つの誘電体の境界条件が設定できないからである．

以上のことから，差分法はこのように誘電率が異なる静電場の計算には適さない．

図 2: 差分法のメッシュとノード．誘電率の異なる部分を横切るので，取り扱いが困難二になる．

2.2 有限要素法

一つの要素で完結すると言うことから，変分法を基礎とする有限要素法の方が，このような問題に適する．変分法とラプラス方程式の関係は，付録 A を見よ． 2.2.1 汎関数の計算

次の汎関数^{1}

J[ϕ] = ∫_{V} \frac{ε}{2} (∇ϕ · ∇ϕ) dV  (12)

の第一変分がゼロの時の，ラプラス方程式になる (付録 A 参照)．この第一変分がゼロになる ϕ が求める— ことが静電場の問題となる．ここで重要なことは，被積分関数が 1 階の微分になっていることである．後で述べることになるが， 1 階の微分だと一つの要素の隣接する 4 つのポテンシャルから計算できる．

汎関数は 3 次元の積分であるが，これ以降，二次元で話を進める．三次元だと図を書くのが大変だし，式も長くなる．また，二次元であろうが三次元であろうが本質的に同じで，三次元への拡張も簡単である．

任意の領域で式 (12) の積分は，図 3 正方形メッシュに分割して近似計算できる．

U[ϕ] = \frac{1}{2} ∑_{i} ∑_{j} ∫_{Ω_{ij}} ε∇ϕ · ∇ϕdV  (13)

四角形要素ごとに積分を行い，すべてを足しあわせることで全体の汎関数 U[ϕ] の値を計算するのである．積分は要素ごとに行うので，誘電率が変化しても要素内で同一になるようにしておけば計算は容易になる．

図 3: 有限要素法のメッシュとノード．おのおの要素ごとに積分を行う．

式 (13) の計算は，四角形要素内での積分を行わなくてはならない．ひとつの四角形要素を図 4 のようにし，それを 4 分割して積分を行う．まず 1 番目の部分の勾配は

∇ϕ = (\frac{ϕ_{i+1 j} − ϕ_{i j}}{h} , \frac{ϕ_{i j+1} − ϕ_{i j}}{h} )  (14)

となる．したがって， 1 番目の積分は，

\frac{1}{2} ∫ ε∇ϕ · ∇ϕdV = \frac{ε_{i+1/2 j+1/2}}{8} [(ϕ_{i+1 j} − ϕ_{i j})^{2} + (ϕ_{i j+1} − ϕ_{i j})^{2}]  (15)

1注意:この汎関数は，静電場のエネルギーになっている．

となる．同じことを 2 番目， 3 番目， 4 番目の領域に対して行い，合計すると i j 番目のエレメントの積分値が計算できる．

∫_{Ω_{ij}} ε∇ϕ · ∇ϕdV = \frac{ε_{i+1/2 j+1/2}}{4} [(ϕ_{i+1 j} − ϕ_{i j})^{2} + (ϕ_{i+1 j+1} − ϕ_{i j+1})^{2}

+(ϕ_{i j+1} − ϕ_{i j})^{2} + (ϕ_{i+1 j+1} − ϕ_{i+1 j})^{2}] (16)

領域全体 Ω の積分は，

U[ϕ] = \frac{1}{2} ∑_{i} ∑_{j} ∫_{Ω_{ij}} ε∇ϕ · ∇ϕdV

= \frac{1}{4} ∑_{i} ∑_{j} ε_{i+1/2 j+1/2} [(ϕ_{i+1 j} − ϕ_{i j})^{2} + (ϕ_{i+1 j+1} − ϕ_{i j+1})^{2}

+(ϕ_{i j+1} − ϕ_{i j})^{2} + (ϕ_{i+1 j+1} − ϕ_{i+1 j})^{2}]  (17)

となる．

図 4: 有限要素法のひとつの i j 番目のエレメントとポテンシャル．要素内で誘電率は一定で， ε_{i j} である．積分は 1〜4 の領域に分けて行う．

2.2.2 第一変分の計算

付録 A で述べたように，式 (17) の個々の ϕ_{i j} を変化させても，汎関数 U[ϕ] の値が変化しないとき，正しいポテンシャル ϕ_{i j} となる．ようするに，境界条件を満たしつつ，静電場のエネルギー U[ϕ] が停留値— ここでは極小値—をとるポテンシャル ϕ を探せということである．

式 (17) 汎関数が停留値になる条件は，

\frac{∂U}{∂ϕ_{1 1}} = 0  \frac{∂U}{∂ϕ_{1 2}} = 0  \frac{∂U}{∂ϕ_{1 3}} = 0  ⋯  \frac{∂U}{∂ϕ_{i j}} = 0  ⋯  (18)

となる．これらの式のうち，もっとも一般的な ∂U/∂ϕ_{i j} を計算する． i = 1, 2, 3, ⋯ ， j = 1, 2, 3, ⋯ と i と j を変化させれば，すべての式を得ることができる．ただし，境界に接する要素は気をつけなくてはならない．

∂U/∂ϕ_{i j} を計算するために，式 (17) の ϕ_{i j} の周りの 4 つの要素に関わる項を書き出すと

U[ϕ] = ⋯

+ \frac{ε_{i−1/2 j−1/2}}{4} [(ϕ_{i j−1} − ϕ_{i−1 j−1})^{2} + (ϕ_{i j} − ϕ_{i−1 j})^{2} + (ϕ_{i−1 j} − ϕ_{i−1 j−1})^{2} + (ϕ_{i j} − ϕ_{i j−1})^{2}]

+ \frac{ε_{i+1/2 j−1/2}}{4} [(ϕ_{i+1 j−1} − ϕ_{i j−1})^{2} + (ϕ_{i+1 j} − ϕ_{i j})^{2} + (ϕ_{i j} − ϕ_{i j−1})^{2} + (ϕ_{i+1 j} − ϕ_{i+1 j−1})^{2}]

+ ⋯

+ \frac{ε_{i−1/2 j+1/2}}{4} [(ϕ_{i j} − ϕ_{i−1 j})^{2} + (ϕ_{i j+1} − ϕ_{i−1 j+1})^{2} + (ϕ_{i−1 j+1} − ϕ_{i−1 j})^{2} + (ϕ_{i j+1} − ϕ_{i j})^{2}]

+ \frac{ε_{i+1/2 j+1/2}}{4} [(ϕ_{i+1 j} − ϕ_{i j})^{2} + (ϕ_{i+1 j+1} − ϕ_{i j+1})^{2} + (ϕ_{i j+1} − ϕ_{i j})^{2} + (ϕ_{i+1 j+1} − ϕ_{i+1 j})^{2}]

+ ⋯  (19)

となる．これは，式 (17) の和の計算の部分を展開しただけであるが，図 5 を見ればこうなることが分かるだろう．この結果を用いると，

\frac{∂U}{∂ϕ_{i j}} = \frac{ε_{i−1/2 j−1/2}}{2} (2ϕ_{i j} − ϕ_{i−1 j} − ϕ_{i j−1}) + \frac{ε_{i+1/2 j−1/2}}{2} (2ϕ_{i j} − ϕ_{i+1 j} − ϕ_{i j−1})

+ \frac{ε_{i−1/2 j+1/2}}{2} (2ϕ_{i j} − ϕ_{i−1 j} − ϕ_{i j+1}) + \frac{ε_{i+1/2 j+1/2}}{2} (2ϕ_{i j} − ϕ_{i+1 j} − ϕ_{i j+1}) = 0 (20)

が得られる．これは，計算領域内部のすべてのポテンシャルについて成り立つ．したがって，式 (20) は， i = 1, 2, 3, ⋯ ， j = 1, 2, 3, ⋯ とすることにより連立方程式となっていることが理解できるだろう．境界条件としてポテンシャルの値が与えれれているところをのぞいて，この連立方程式を解けばよい．

図 5: 汎関数の積分計算のうち，ポテンシャル ϕ_{i, j} が関係する要素．

領域の境界でポテンシャルが与えられていない場合，その場所では電場が境界と垂直になる^{2}．これは，式 (25) から保証される．

連立方程式 (20) を見ると，ほとんど差分の式と同じである．もし，誘電率が一定とすると差分の式と全く同一になる．

2.3 有限積分法

もちろん，有限積分法も適用できる．時間がないので，その解説は行わない．

2いわゆるノイマン条件

付録 A 変分法

計算領域 Ω でのポテンシャルを ϕ とする．このポテンシャルをパラメーターとした汎関数 U[ϕ]

U[ϕ] = ∫_{V} \frac{ε}{2} (∇ϕ · ∇ϕ) dV  (21)

を考える．これから，この汎関数が停留値をとるとき，そのときの ϕ はラプラス方程式を満足することを示す．ラプラス方程式を満足するので，真の解となる．

汎関数が停留値になる条件は，その第一変分がゼロである．式 (21) の凡関数の第一変分 δU は

δU = \{ lim_{α→0} \frac{U[ϕ + αv] − U[ϕ]}{α}\} α  (22)

である．ここで α は任意の実数， v は領域 Ω の境界ではゼロとなる任意の関数である．したがって，汎関数の式 (21) の第一変分は，

δU = \{ lim_{α→0} \frac{∫_{Ω} \frac{ε}{2}∇(ϕ + αv) · ∇(ϕ + αv)dV − ∫_{Ω} \frac{ε}{2}∇ϕ · ∇ϕdV}{α} \} α
= α ∫_{Ω} ε∇ϕ · ∇vdV  (23)

となる．ここで，次の微分

∇ · (εv∇ϕ) = ε∇v · ∇ϕ + v∇ · (ε∇ϕ)  (24)

を使うと，汎関数式 (23) は

δU = α ∫_{Ω} ∇ · (εv∇ϕ)dV − α ∫_{Ω} v∇ · (ε∇ϕ)dV
= α ∫_{Ω} εv∇ϕ · ndS − α ∫_{Ω} v∇ · (ε∇ϕ)dV  (25)

となる．右辺第一式の積分はゼロとなる．なぜならば，関数 v は考えてる領域 Ω の境界でゼロとなる関数であるからである．したがって，

δU = −α ∫_{Ω} v∇ · (ε∇ϕ)dV  (26)

である．汎関数の第一変分 δU がゼロとなるためには，

∇ · (ε∇ϕ) = 0  (27)

とならなくてはならない．関数 v は任意の関数であるため， ∇ · (ε∇ϕ) がゼロの場合のみ積分の値がゼロになる． ∇ · (ε∇ϕ) がゼロでなくても特定の v で積分がゼロになることはある．しかし，任意の v となると話は別である．

ここでの結論は，「式 (21) の汎関数の第一変分をゼロにすることと，式 (27) の微分方程式と解くことは等価」である．ようするに，式 (27) の微分方程式の解は，式 (21) の積分が停留値—しばしば極小値—になるような関数 ϕ と同一である．

式 (21) の積分は，静電場のエネルギーになっている．境界条件を満たしつつ，このエネルギーが最低になる状態が実際の状態である．ここでは，この積分が極小になる ϕ を計算することになる．

この静電場も最小エネルギーの問題になっている．しばしば，物理的に安定な状態はエネルギーが最低になる．この静電場の問題もその例の一つである．