離散化した場合の微分 (回転と勾配)

山本昌志^{∗}

2005 年 2 月 9 日

1 はじめに

マクスウェルの方程式を離散化して、有限要素法により、電磁場を計算しようとしている。実際には、電磁場そのものではなく、ポテンシャルを計算することが多い。電磁場を知りたければ、それを微分すれば良い。有限要素法の計算では、離散化されたポテンシャルが得られるが、その微分について少し考えることにする。ここでは、離散化された量の微分の計算方法を示す。

今のところ我々は、軸対称形状の以下の問題を計算している。

• 静電場の問題

• 静磁場の問題

• 共振空洞の問題

これらの電場や磁場を求める場合の微分の方法について示す。ただし、ここでは、軸対称問題視か扱わないので、円柱座標系を用いることにする。円柱座標系の詳細を知りたければ、「座標系と微分演算子^{1}」を見よ。

2 回転

回転から電磁場を求めるのは、静磁場と共振空洞の問題に現れる。即ち、

B = ∇ × A  (静磁場)  (1)
E = \frac{i}{ε_{0}ω} ∇ × H  (共振空洞)  (2)

である。微分 (回転) と言う意味では、全く同じ形をしている。そのため、ここでは静磁場の問題のみについて述べる。共振空洞の場合は同じように考えて欲しい。

∗独立行政法人秋田工業高等専門学校 電気情報工学科

1http://www.akita-nct.jp/ yamamoto/study/electromagnetics/coodinate transform/html/index.html

2.1 要素内

軸対称構造、実際には軸対称電磁石であるが、その場合の有限要素解析では、ベクトルポテンシャルを計算するのが一般的である。軸対称なので、ベクトルポテンシャルは θ 方向成分のみである。従って、このベクトルポテンシャルの回転から求められる磁場は、

H_{r} = −\frac{∂A_{θ}}{∂z}  (3)
H_{z} = \frac{1}{r} \frac{∂}{∂r} (rA_{θ})  (4)

となる。

有限要素法では、図 1 に示すように、三角形要素の頂点でベクトルポテンシャル A_{θ} の値を計算する。これら 3 つのベクトルポテンシャルから、要素内の磁場 H_{r} と H_{z} を求めなくてはならない。まず初めに、要素内の任意の位置でのベクトルポテンシャルを計算する式を導く。今のところ、我々の計算は 1 次近似なので、要素内のベクトルポテンシャルは

A_{θ}(z, r) = α + βz + γr  (5)

となる。

この要素内のベクトルポテンシャルの値を用いて、領域内の磁場を計算することにする。これは簡単で、式 (5) を用いて、式 (3) と (4) を評価すれば良い。これらから

H_{r} = −β  (6)
H_{z} = \frac{α + βz}{r} + 2γ  (7)

を導くことができる。ここで、 H_{r} の式に座標 z が含まれるのが奇異に感じるかもしれない。これは、要素内でベクトルポテンシャル A_{θ} が一定でも、回転が生じることを考えると理解できる。

残りの問題は、ベクトルポテンシャルを 1 次近似した係数 α と β、 γ を求めることである。これが分かれば、要素内の任意の位置での磁場が計算できる。三角形要素の頂点の座標とベクトルポテンシャルの有限要素法の計算結果によりが分かっているので、これら 3 つの係数を決定することは可能である。これら 3 点では

[\begin{matrix}1 & x_{i} & y_{i} \\ 1 & x_{j} & y_{j} \\ 1 & x_{k} & y_{k} \end{matrix}] [\begin{matrix}α \\ β \\ γ \end{matrix}] = [\begin{matrix}A_{θ i} \\ A_{θ j} \\ A_{θ k} \end{matrix}]  (8)

の関係がある。ここで、 (x_{i}, y_{i}) と (x_{j}, y_{j})、 (x_{k}, y_{k}) は頂点の座標、 A_{θ i} と A_{θ j}、 A_{θ k} はその位置のベクトルポテンシャルである。この連立方程式から、 1 次近似の係数は

[\begin{matrix}α \\ β \\ γ \end{matrix}] = [\begin{matrix}1 & x_{i} & y_{i} \\ 1 & x_{j} & y_{j} \\ 1 & x_{k} & y_{k} \end{matrix}]^{−1} [\begin{matrix}A_{θ i} \\ A_{θ j} \\ A_{θ k} \end{matrix}]  (9)

と計算できる。これで、要素内の磁場の値は計算できる。図 1: 三角形要素とベクトルポテンシャル

3

2.2 内部要素の頂点

要素内の磁場は、式 (6) や (7) を用いて計算できる。しかし、エレメントの境界での値を考えると問題が生じる。例えば、図 2 に示すように要素 A,B,C,D,E,F が共有する頂点の磁場の値をどうするかである。今のところ、このような場合は平均とするのが良いであろう。

図 2: 隣接要素の様子

2.3 軸上

式 (7) を見ると、軸上 (r = 0) では問題が起きそうである。ゼロで割るため発散する可能性がある。問題となるのは、図 3 に示す要素うち、 A や C のように要素の一辺が軸となる場合と B のように 1 つの頂点が軸上にある場合である。

まずは、要素の一辺が軸となる場合について考える。この要素の軸上のベクトルポテンシャルの値は、物理的にゼロ (A_{θ} = 0) にならなくてはならない。そのため、要素内のベクトルポテンシャルを示す式 (5) は、

0 = α + βz  z_{1} ≦ z ≦ z_{2}  (10)

となる。図 3 の要素 A の場合、軸となる要素の一辺 (z_{1} ≦ z ≦ z_{2}) で、これが成立しなくてはならないので、

α = 0  (11)
β = 0  (12)

となる必要がある。要するに、要素の一辺が軸上にある場合、その α と β はゼロになり、ベクトルポテンシャルは r 座標に比例しているのである。そのため、式 (7) が発散する事はなく、問題は生じないのである。まとめると、この場合の軸上の磁場は、

H_{r} = 0  (13)
H_{z} = 2γ  (14)

となる。軸上では H_{z} のみで、これは物理的に正しい。非常にめでたい。

次に、図 3 の要素 B のように、頂点の一つが軸上にある場合を考える。この軸上の頂点でのベクトルポテンシャルはゼロなので、

0 = α + βz_{2}  (15)

である。この式は頂点での一点 (r = 0, z = z_{2}) でしか成り立たない。そのため、一辺が軸上にあるときのような式 (11) や (12) に対応するものは無い。従って、この要素の軸上の頂点での磁場は

H_{r} = −2β  (16)
H_{z} = 2γ  (17)

となる。 H_{z} が発散しないのは良いが、 H_{r} がゼロにならないのはまずい。

以上のことから、図 3 のように軸上のノードの磁場は、図中の要素 A や B のように一辺が軸上にある磁場の平均とするのが良いだろう。 H_{z} の発散の問題も生じないし、 H_{r} が有限になることもない。

図 3: 軸付近の要素の様子

2.4 境界

軸上の磁場の計算で示したように、頂点の一つが軸上にある場合、物理的な条件に合致しない場合がある。先の場合は、軸上にある一つの頂点での H_{r} がゼロとならないのである。一方、一辺が軸上にある場合は、物理的な条件に完全に合致している。そのような理由から、軸上にあるノードでの磁場の値は、それを共有するものうち、軸上に辺がある要素の平均を取ることにした。

これと同じことが、軸以外の境界でも生じる。特に、境界が直線の場合は、全く同じ状況になるであろう。このことから、ここでも軸と同じように考えるのが良いだろう。即ち、図 4 の様な境界があるとき、要素 A,B,C が共有するノードの磁場の値は、要素 A と B の平均とする。要素 B は考えないものとする。

図 4: 境界付近の要素の様子

3 勾配

勾配から電磁場を求めるのは、静電場の問題に現れる。即ち、

E = −∇ϕ  (18)

である。ここではこの勾配の計算方法について述べる。

3.1 要素内

静電場の問題はスカラーポテンシャルを計算するのが一般的である。それがもっとも計算が簡単であるからである。我々は、軸対称構造の静電場の問題を計算している。軸対称なので、ポテンシャル ϕ は (r, z) の関数で、 θ 方向には依存しない。従って、このポテンシャルの勾配から求められる電場は、

E_{r} = −\frac{∂ϕ}{∂r}  (19)
E_{z} = −\frac{∂ϕ}{∂z}  (20)

となる。

本図書の回転で示したのと同様に、有限要素法では図 1 で示される三角形要素の頂点でスカラーポテンシャル ϕ の値を計算する。これら 3 つのスカラーポテンシャルから、要素内の電場 E_{r} と E_{z} を求めなくてはならない。要素内の任意の位置でのスカラーポテンシャルをを計算する式を導く。今のところ、我々の計算は 1 次近似なので、要素内のスカラーポテンシャルは

ϕ(z, r) = α + βz + γr  (21)

となる。ポテンシャルの値が分かれば、電場は計算できる。式 (21) を式 (19) や (20) に代入すると、

E_{r} = −γ  (22)
E_{z} = −β  (23)

となる。スカラーポテンシャルの値を示す 1 次近似の係数 γ や β を計算すれば、電場が分かるのである。非常に単純な式である。

電場を計算するための係数 α や γ、 β の計算は、回転の場合と同様にすればよく、

[\begin{matrix}α \\ β \\ γ \end{matrix}] = [\begin{matrix}1 & x_{i} & y_{i} \\ 1 & x_{j} & y_{j} \\ 1 & x_{k} & y_{k} \end{matrix}]^{−1} [\begin{matrix}ϕ_{i} \\ ϕ_{j} \\ ϕ_{k} \end{matrix}]  (24)

となる。もちろん、 (x_{i}, y_{i}) と (x_{j}, y_{j})、 (x_{k}, y_{k}) は頂点の座標、 ϕ_{i} と ϕ_{j}、 ϕ_{k} はその位置のスカラーポテンシャルである。これで、要素内の電場の値は計算できる。

3.2 内部要素の頂点

内部要素で共有された頂点での電場の値を求めるときには、回転の時と同様に問題が生じる。図 2 のように要素 A,B,C,D,E,F が、頂点 i を共有し、その点での電場の値をどうするかである。 A〜F のそれぞれの要素で電場の値が微妙に異なることが問題である。その差は小さいので、とりあえず平均値とするのが良いであろう。

3.3 軸上

式 (22) や (23) を見ても分かるように、回転の時のように値が発散するようなことはない。しかし、境界条件について、少し考えてみるのも良かろう。軸対称問題の静電場の場合、その軸上では、自然境界条件

(∇ϕ) · n = 0  (25)

を満たす必要がある。ここで、 n は軸に垂直で外側に向かっている単位ベクトルである。これは、 −r 方向を向いている。 ∇ϕ の −r 方向成分は、 −∂ϕ/∂r = γ である。従って

γ = 0  (26)

となる必要がある。

軸上の要素、例えば図 3 の要素 A, B, C の γ がそのようになっているかは不明である。計算精度にも依存するであろう。従って、ここでは変な小細工をしないで、軸上でも E_{r} が出てきても、それは計算誤差として考えるのが良いであろう。このことから、我々は回転の計算同様に、図 3 の要素 A, B, C が共有するノードでの電場 (E_{r} と E_{z}) は、要素 A と C の平均とする。要素 B のように軸上の 1 つのノードしかない要素は、平均の計算に入れないこととする。

3.4 境界

この場合も回転と同様に考える。即ち、図 4 の様な境界があるとき、要素 A,B,C が共有するノードの電場の値は、要素 A と B の平均とする。要素 B は考えないものとする。