Astroの光線のサムネイル。

pubDate: 2024-03-18

author: sakakibara

有限要素法

電磁場解析

数値解析

電磁界解析

マクスウェル方程式を手で解く人はいない。
解析的に解けないケースがほとんどだからだ。 実際に電磁界解析を行おうと思うとコンピュータによる数値解析が必須になってくる。 だが、電磁界解析を行うソフトは非常に高価で数100万~数1000万程度の金かかる。 それにカスタマイズ性も低く、使いこなすのに練度が必要となる。

そこで、自前の電磁界解析プログラムを作成して、様々な実験を行う。

以下ではどのようにしてマクスウェル方程式を数値解析に落とし込むか。 より具体的に言えばどのようにして微分方程式を行列計算に落とし込むかについて説明する。

重み付き残差法

内積

内積という概念が広いことを知っている人は多いだろう。

⟨⋅⟩:V×V→R\langle\cdot \rangle: V\times V \rightarrow \mathbb{R} \\

a,b,c∈Va, b, c \in V, α∈R\alpha \in \mathbb{R}について以下が成り立つなら演算子⟨⋅⟩\langle\cdot\rangleを内積という。

(正値性) ⟨a⋅a⟩≥0 ∧ ⟨a⋅a⟩=0  ⟹  a=0(対称性) ⟨a⋅b⟩=⟨b⋅a⟩(多重線形性) α⟨a⋅b⟩=⟨αa⋅b⟩ ⟨a+b⋅c⟩=⟨a⋅c⟩+⟨b⋅c⟩\begin{aligned} \text{(正値性)}&\ \langle a\cdot a \rangle \ge 0\ \land\ \langle a\cdot a \rangle = 0 \implies a = 0\\ \text{(対称性)}&\ \langle a\cdot b \rangle = \langle b \cdot a \rangle \\ \text{(多重線形性)}&\ \alpha \langle a \cdot b \rangle = \langle \alpha a \cdot b \rangle \\ &\ \langle a + b \cdot c \rangle = \langle a \cdot c \rangle + \langle b \cdot c \rangle \end{aligned}

複素数を前提とした定義もあるが、あまり使用しないので今回は以上を満たす演算⟨⋅⟩\langle \cdot \rangleを内積とする。

内積という概念は距離の概念を拡張したものと取られられている。 これは関数としての距離が内積の定義に含まれるからである。

残差を求める方法

ベクトル空間VV、f:V→Vf : V\rarr Vとする。
a,b∈V\bm{a}, \bm{b} \in Vにたいして

f(a)=bf(\bm{a}) = \bm{b}

以上のような関数を考える。ここでf,bf, \bm{b}を既知としてa\bm{a}を求めたいがどのようにすればよいだろうか。
最初に考得られる手段が最小二乗法である。
a′\bm{a}'をa\bm{a}の近似解とする。以下で定義するR\bm{R}を残差という。

R=b−f(a′)\bm{R} = \bm{b} - f(\bm{a}')

最小二乗法は残差の二乗和(内積)を最小にするような近似解a′\bm{a}'をa\bm{a}とする手法である。

a=min⁡a′(⟨R(a′)⋅R(a′)⟩)\bm{a} = \min_{\bm{a'}} \big(\langle\bm{R}(\bm{a}')\cdot\bm{R}(\bm{a}')\rangle\big)

残差の二乗を最小にするベクトルが解であるというのは非常に直感的だが、この他にこれに近い考えとして重み付き残差法というのものがある。

0=⟨w⋅R(a′)⟩0 = \langle\bm{w}\cdot\bm{R}(\bm{a}')\rangle

ここでw\bm{w}は重み関数と呼ばれる関数で、この種類によって様々なバリエーションが生まれる。

ガラーキン法

ベクトル空間VVが基底⟨N1,N2,…,NE⟩\langle \bm{N}_ 1, \bm{N}_ 2, \ldots, \bm{N}_ E \rangleによって生成されているとする。 a∈V\bm{a}\in Vについて以下が成り立つ。

a=∑aeNe\bm{a} = \sum a_e \bm{N}_ e

そこで、重み付き残差法におけるw\bm{w}をNe\bm{N}_ eとすると、

0=⟨Ne⋅R(a′)⟩0 = \langle\bm{N}_ e\cdot\bm{R}(\bm{a}')\rangle

が得られる。
Ne\bm{N}_ eがEE個あるため方程式がEE個得られる。 これによりEE個のaea_eを変数とする方程式が得られるため、行列で計算できる。 ガラーキン法の特徴として残差が基底により生成する空間の直交空間となることである。

有限要素法

有限要素法とは対象(これをMMとする。)を複数の要素(または有限要素)と呼ばれる領域に分割して、各領域で物理量を補間することにより実際の現象をシミュレーションする手法のことだ。 分割された要素は22次元であれば三角形や四角形であったり、33次元であれば四面体、六面体、プリズム、などに分割される。基本的に分割が細ければ細かいほどシミュレーションとしての精度は向上するが、計算コストが増大する。
また、要素内における物理量の補間に関しても1次補間、2次補間、…といった分類やラグランジュ補間、ルジャンドル補間などの補間手法による分類などががある。

今回は最も単純な四面体要素、1次ラグランジュ補間について解説する。

A-ϕ\phi法

よく知られたようにマクスウェル方程式は以下の４つの式で表される。

divD=ρdivB=0rotE=−∂B∂trotH=J+∂D∂t\begin{align} \mathrm{div}\bm{D} &= \rho \\ \mathrm{div}\bm{B} &= 0 \\ \mathrm{rot}\bm{E} &= -\frac{\partial \bm{B}}{\partial t} \\ \mathrm{rot}\bm{H} &= \bm{J}+\frac{\partial \bm{D}}{\partial t} \end{align}

ただし、適当なテンソルϵ,μ,σ\epsilon, \mu, \sigmaを用いて

D=ϵEB=μHJ=σE+Jsource\bm{D} = \epsilon \bm{E} \\ \bm{B} = \mu \bm{H} \\ \bm{J} = \sigma \bm{E} + \bm{J}_ {source}

のような関係がある。

磁束密度B\bm{B}はベクトル・ポテンシャルにより表現することもできる。 ベクトル・ポテンシャルをA\bm{A}とすると。

B=rot(A+gradχ)\bm{B} = \mathrm{rot}(\bm{A} + \mathrm{grad}\chi)

のように表現できる。実質的にB\bm{B}の値を決定しているのはA\bm{A}であるが、 rot\mathrm{rot}の演算はgrad\mathrm{grad}を00にするためにgradχ\mathrm{grad}\chiの分だけ自由度が残る。 これをゲージ自由度という。 これを33つめの式に代入すると

rotE=−∂B∂t=−∂∂trot(A+gradχ)=−rot(∂A∂t+∂∂tχ)\begin{align} \mathrm{rot}\bm{E} &= -\frac{\partial \bm{B}}{\partial t} \\ &= -\frac{\partial}{\partial t}\mathrm{rot}(\bm{A} + \mathrm{grad}\chi) \\ &= -\mathrm{rot}(\frac{\partial \bm{A}}{\partial t} + \frac{\partial}{\partial t}\chi) \end{align}

より

rot(E+∂A∂t+∂∂tχ)=0\mathrm{rot}(\bm{E} + \frac{\partial \bm{A}}{\partial t} + \frac{\partial}{\partial t}\chi) = \bm{0}

これにより(ヘルムホルツの原理)、

E+∂A∂t+∂∂tχ=gradχ′\bm{E} + \frac{\partial \bm{A}}{\partial t} + \frac{\partial}{\partial t}\chi = \mathrm{grad}\chi' E=−∂A∂t−∂∂tχ+gradχ′\bm{E} = -\frac{\partial \bm{A}}{\partial t} - \frac{\partial}{\partial t}\chi + \mathrm{grad}\chi'

ここでもχ,χ′\chi, \chi'の分だけ不定となる。 この不定を決定することによりさまざまなパターンのベクトル・ポテンシャルが定義できる。 不定の値を様々にチューニングしてベクトル・ポテンシャル間を変換することをゲージ変換という。 逆にマクスウェル方程式はゲージ変換の元で普遍な理論という見方もできる。 (実用上、計算がしやすくなる以上の意味は見いだせない。)

以降の有限要素法ではχ=0\chi=0, χ′=∂ϕ∂t\chi' = \frac{\partial \phi}{\partial t}とする。つまり

B=rotAE=−∂A∂t+grad∂ϕ∂t\begin{aligned} \bm{B} &= \mathrm{rot}\bm{A} \\ \bm{E} &= -\frac{\partial \bm{A}}{\partial t} + \mathrm{grad}\frac{\partial \phi}{\partial t} \end{aligned}

を選択する。 したがって。

J=σE+Jsource=−σ∂A∂t−gradσ∂ϕ∂t+Jsource\begin{aligned} \bm{J} &= \sigma\bm{E} + \bm{J}_ {source} \\ &= -\sigma\frac{\partial \bm{A}}{\partial t} - \mathrm{grad}\sigma\frac{\partial \phi}{\partial t} + \bm{J}_ {source} \end{aligned}

少しだけ注意すべきはE\bm{E}の第二項が∂ϕ∂t\frac{\partial \phi}{\partial t}となっていることである。 多くの参考資料ではこの代わりにE=−∂A∂t+ϕ\bm{E} = -\frac{\partial \bm{A}}{\partial t} + \bm{\phi}となっていることがおおいだろう。

さて、これまでのA,∂ϕ∂t\bm{A}, \frac{\partial \phi}{\partial t}の導入から式33を満たすことは自明である。 式44について考える。

rotH=J+∂D∂trotμ−1B=J+ϵ∂E∂trotμ−1rotA=J−ϵ∂2A∂t2+ϵgrad∂2ϕ∂t2\begin{aligned} \mathrm{rot}\bm{H} &= \bm{J}+\frac{\partial \bm{D}}{\partial t} \\ \mathrm{rot}\mu^{-1}\bm{B} &= \bm{J}+\epsilon\frac{\partial \bm{E}}{\partial t} \\ \mathrm{rot}\mu^{-1}\mathrm{rot}\bm{A} &= \bm{J}-\epsilon\frac{\partial^2 \bm{A}}{\partial t^2} + \epsilon\mathrm{grad}\frac{\partial^2 \phi}{\partial t^2} \\ \end{aligned}

これ以上の式変形は不要に複雑になるのでここで留めておこう。

有限要素法では右辺第二項以降を無視することが多い。 これは準静的過程(quasi-static)と呼ばれる。英語名が不思議なので検索に意外と引っかからない。 第二項を無視できる条件は周波数と誘電率と伝導率が関わるが、大体無視できるのでよい。

rotμ−1rotA=Jrotμ−1rotA=−σ∂A∂t−σgrad∂ϕ∂t+Jsource\begin{aligned} \mathrm{rot}\mu^{-1}\mathrm{rot}\bm{A} &= \bm{J} \\ \mathrm{rot}\mu^{-1}\mathrm{rot}\bm{A} &= -\sigma\frac{\partial \bm{A}}{\partial t} - \sigma\mathrm{grad}\frac{\partial \phi}{\partial t} + \bm{J}_ {source} \\ \end{aligned}

少し整理して

rotμ−1rotA+σ∂A∂t+σgrad∂ϕ∂t=Jsource\begin{equation} \mathrm{rot}\mu^{-1}\mathrm{rot}\bm{A} + \sigma\frac{\partial \bm{A}}{\partial t} + \sigma\mathrm{grad}\frac{\partial \phi}{\partial t} = \bm{J}_ {source} \tag{eq:$A-\phi$} \\ \end{equation}

有限要素法ではこの式(eq:A−ϕ)(eq:A-\phi)を解く。

補間と形状関数

位置次元R\mathbb{R}上で定義された実数値関数f(x)f(x)について, f(x1),f(x2),…,f(xn)f(x_1), f(x_2), \ldots, f(x_n)が与えれているとする。この元の関数f(x)f(x)はどのような関数だろうか？

単純な考えとして、実数値関数f(x)f(x)を多項式pn−1(x)p_{n-1}(x)で近似する方法だ。

pn−1(x)=a0+a1x+a2x2+a3x3⋯+an−1xn−1≃f(x)\begin{aligned} p_{n-1}(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 \cdots + a_{n-1} x^{n-1} \\ &\simeq f(x) \end{aligned}

このような関数を補間する多項式のことを補間多項式と呼ぶ。 のように表現して連立方程式を解く方法が考えられる。

(f(x1)f(x2)⋮f(xn))=(1x1x12⋯x1n−11x2x22⋯x2n−1⋮1xn−1xn−12⋯xn−1n−1)(a0a1⋮an−1)\begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{pmatrix} = \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ &&&\vdots& \\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{n-1}^{n-1} \\ \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{pmatrix}

この真ん中の行列はヴァンデルモンドの行列として行列式が非常にきれいに求められることが知られている。(可逆である。) この解法に見られる特徴として、f(xi)=p(xi)f(x_i)=p(x_i)となっている点である。これは補間という考えで非常に重要な性質である。

実際にこれは最もシンプルな方法であり、計算もうまくことが知られている。 ただし、これは計算コストが高いことが知られている。(LU分解ができない)

もっとうまいやり方がある。

それは以下のような補間多項式を用いる方法だ。

p2(x)=f(x1)L1(x)+f(x2)L2(x)+f(x3)L3(x)p_{2}(x) = f(x_1)L_1(x) + f(x_2)L_2(x) + f(x_3)L_3(x)

ただし、LiL_iは22次の多項式であり以下を満たす。

Li(xj)={1if (i=j)0if (i≠j)L_i(x_j) = \begin{cases} 1 & \text{if } (i = j) \\ 0 & \text{if } (i \neq j) \end{cases}

補間多項式は33つの項で表されているが、実際には22次の多項式となる。 この関数ϕi\phi_iを用いる補間をラグランジュ補間と呼ぶ。
これをヴァンデルモンドを用いる方法で再び示してみよう。

(p(x1)p(x2)⋮p(xn))=(L1(x1)0⋯00L0(x0)⋯0⋮⋮00⋯Ln(xn))(f(x1)f(x2)⋮f(xn))\begin{pmatrix} p(x_1) \\ p(x_2) \\ \vdots \\ p(x_n) \end{pmatrix} = \begin{pmatrix} L_1(x_1) & 0 & \cdots & 0 \\ 0 & L_0(x_0) & \cdots & 0 \\ \vdots&&&\vdots&\\ 0 & 0 & \cdots & L_{n}(x_n)\\ \end{pmatrix} \begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{pmatrix}

このようにヴァンデルモンドを使用していた部分が単位行列になる。
このラグランジュ補間はnn個の標本点に対してn−1n-1次の多項式を構成する。 最終的に得られる補間多項式は形式的には(展開すれば)、上で挙げたヴァンデルモンドの方程式と変わらない。 Li(xj)L_i(x_j)のように標本点で11となるようなn−1n-1次の多項式を基底関数に用いた補間多項式は一意に定まる。 つまり、ヴァンデルモンドの行列を計算して導出しても、ラグランジュ補間を用いて導出した補間多項式しても、結果的に得られる補間多項式は同じである。
今回の場合は、これらは具体的に以下のように書き下せる。

L1(x)=(x−x2)(x−x3)(x1−x2)(x1−x3)L2(x)=(x−x1)(x−x3)(x2−x1)(x2−x3)L3(x)=(x−x1)(x−x2)(x3−x1)(x3−x2)\begin{aligned} L_1(x) = \frac{(x - x_2)(x - x_3)}{(x_1 - x_2)(x_1 - x_3)} \\ L_2(x) = \frac{(x - x_1)(x - x_3)}{(x_2 - x_1)(x_2 - x_3)} \\ L_3(x) = \frac{(x - x_1)(x - x_2)}{(x_3 - x_1)(x_3 - x_2)} \end{aligned}

また、これらの関数は互いに直交することがわかる。(部分積分をしてみればすぐに気づく。どのみちxix_iを代入すれば00になるため)

まとめる。
一般にn−1n-1次の多項式は

pn−1(x)=f(x1)L1(x)+f(x2)L2(x)+⋯+f(xn)Ln(x)p_{n-1}(x) = f(x_1)L_1(x) + f(x_2)L_2(x) + \cdots + f(x_{n})L_{n}(x)

ただし、Li(xj)L_i(x_j)は以下を満たすn−1n-1次多項式である。

Li(xj)={1if (i=j)0if (i≠j)L_i(x_j) = \begin{cases} 1 & \text{if } (i = j) \\ 0 & \text{if } (i \neq j) \end{cases}

標本点を多く取れば取るほど精度が高くなると考えられる。しかし、ラグランジュ補間は急峻な関数を近似する際、標本点を多く取ると多項式補間が振動する現象(ルンゲ現象)を引き起こすことが知られている。

(なんと、この関数は1779年に考え出されたものらしい。しかも発案者はラグランジュではなく、エドワード・ワーリングという人らしい。)

要素内補間

対象MMのある33次元四面体要素をmmとする。 mm内の物理量ϕ\phiを補間する方法について考える。予め名前をつけておくが、要素内の点における物理量はたいていある関数の線形和で表される。その関数を形状関数とか基底関数などと呼ばれる。 mmは44つの頂点と66つの辺と44つの面を持つ。 44つの頂点にϕ(x0),ϕ(x1),ϕ(x2),ϕ(x3)\phi(\bm{x}_ 0), \phi(\bm{x}_ 1), \phi(\bm{x}_ 2), \phi(\bm{x}_ 3)が割てられているとする。この定数をそれぞれϕ0,ϕ1,ϕ2,ϕ3\phi_0, \phi_1, \phi_2, \phi_3とする。

四面体全体の体積をVVとし、頂点iiを除く頂点と要素mm内の点x\bm{x}で構成される四面体の体積Vi(x)V_i(\bm{x})を考える。

Li(x)=Vi(x)VL_i(\bm{x}) = \frac{V_i(\bm{x})}{V}

とする。 このようにLi(x)L_i(\bm{x})を定義することで

Li(xj)={1if (xi=xj)0if (xi≠xj)L_i(\bm{x}_ j) = \begin{cases} 1 & \text{if } (\bm{x}_ i = \bm{x}_ j) \\ 0 & \text{if } (\bm{x}_ i \neq \bm{x}_ j) \end{cases}

を実現できる。

ϕ(x)=ϕ0L0(x)+ϕ1L1(x)+ϕ2L2(x)+ϕ3L3(x)\phi(\bm{x}) = \phi_0L_0(\bm{x}) + \phi_1L_1(\bm{x}) + \phi_2L_2(\bm{x}) + \phi_3L_3(\bm{x})

を考える。

ここで44つの標本点(x0,ϕ0),(x1,ϕ1),(x2,ϕ2),(x3,ϕ3)(\bm{x_0}, \phi_0), (\bm{x_1}, \phi_1), (\bm{x_2}, \phi_2), (\bm{x_3}, \phi_3)があるなら3次の多項式が得られるのではないか？という疑問が浮上する。

実際に、最初の考えに立ち返ってみる。

ϕ(x)=a0+ax+by+cz\phi(\bm{x}) = a_0 + ax + by + cz

ここでa0,a,b,ca_0, a, b, cの44つを変数とするために44つの式が必要となる。 そしてそれは44つの標本点によって構成できる。 そう。つまり単純にnn個の標本点があればn−1n-1次補完多項式を定義できるというわけではない。 多次元の空間を考える場合には少し複雑になる。

以上によって11次式で表されるLi(x)L_i(\bm{x})を用いてϕ(x)\phi(\bm{x})をラグランジュ補間することができる。

頂点iiに標本点をもつ補間の基底関数を節点形状関数ωi\omega_iという。 この他にも辺eeに標本点をもつ補間の基底関数を辺形状関数Ne\bm{N}_ e、面ffに標本点をもつ補間の基底関数を面形状関数MfM_fという。 そしてこれら形状関数は以下を要請する。

ωi(xj)=δij∫eNedle′=δee′∬fMfdSf′=δff′\begin{aligned} \omega_ i({\bm{x}_ j}) &= \delta_{ij} \\ \int_e\bm{N}_ e\mathrm{d}\bm{l}_ {e'} &= \delta_{ee'} \\ \iint_f M_f \mathrm{d}\bm{S}_ {f'} &= \delta_{ff'} \end{aligned}

なぜ辺形状関数などを使うのかというと連続に関する制限を締め付けないためにある。 頂点形状関数のみだと要素の節点での値が連続であることを要請する。 しかし辺形状関数だとその辺の対角にある面(その辺を含まない面)でのみ接線方向で連続となる。 この性質は媒質の境界面で物性が階段的に変化するような対象をシミュレーションする際に有効となる。 たとえば、空気と鉄を有限要素にりシミュレーションする際、空気の要素と鉄の要素が面をまたいで要素分割されるのであれば、面の垂直方向(面の法線方向)での連続性は考慮しなくて済む。

(さすがに文字だけではイメージできないので図を載せたい。)

面形状関数についてはよく知らない。

ともかく以上のような理由から電磁ベクトル・ポテンシャルの解析には辺形状関数を用いることが多い。

そこで要素mm内におけるベクトルポテンシャルA\bm{A}を次のように辺形状関数Ne\bm{N}_ eを用いて補間したい。

A(x)=∑e∈EdgeAeNe(x)\bm{A}(\bm{x}) = \sum_{e\in \mathrm{Edge}}A_e\bm{N}_ e(\bm{x})

ベクトル・ポテンシャルはその法線方向(垂直方向)には一般には不連続性を要求するからである。 (これは異なる透磁率をもつ物体において磁束密度B\bm{B}が境界面に垂直方向で連続であり、接線方向では一般に不連続であることからくる。B=rotA\bm{B} = \mathrm{rot}\bm{A}として表されるのでベクトルA\bm{A}は境界面の垂直方向には不連続であり、境界面の接線方向には連続である。(ストークスの定理より))

電磁ポテンシャル(まぁ、電位)は節点形状関数LiL_iを用いて補間する。

ϕ(x)=∑n∈NodeϕiLi(x)\phi(\bm{x}) = \sum_{n \in {Node}}\phi_i L_i(\bm{x})

ここで何食わぬ顔でωi\omega_iをテクニカルに微分しよう。

∇ωi=∇ωi∑ωj−ωi∇∑ωj=∑j(ωj∇ωi−ωi∇ωj)\begin{aligned} \nabla \omega_i &= \nabla \omega_i \sum \omega_j - \omega_i \nabla \sum \omega_j \\ &= \sum_j (\omega_j \nabla \omega_i - \omega_i \nabla \omega_j) \end{aligned}

ここで勾配・回転・発散の資料を参考にしてほしいのだが、勾配行列GGを用いると

∇ωi=∑GkiNk\nabla \omega_i = \sum G_{ki} \bm{N}_ k

が得られる。 そこで頂点iiから頂点jjへ向かう辺kkの辺形状関数Nk\bm{N}_ kを 次のように定義する。

Nk=ωi∇ωj−ωj∇ωi\bm{N}_ k = \omega_i \nabla \omega_j - \omega_j \nabla \omega_i

このように定義すると

∇ωi=∇ωi∑ωj−ωi∇∑ωj=∑j−(ωi∇ωj−ωj∇ωi)=∑kGkiNk\begin{aligned} \nabla \omega_i &= \nabla \omega_i \sum \omega_j - \omega_i \nabla \sum \omega_j \\ &= \sum_j -(\omega_i \nabla \omega_j - \omega_j \nabla \omega_i) \\ &= \sum_k G_{ki}\bm{N}_ k \end{aligned}

第二項は辺kkが頂点iiから出ていって、頂点jjへ入っているのでこのようになる。
そしてこれは勾配行列と頂点形状関数の関係を充たす。

(辺形状関数の導出がテクニカルすぎるって？ 俺もそーおもう)

辺形状関数のrot\mathrm{rot}については辺形状関数を以上のように定義したことにより非常に計算しやすくなる。

rotNk=rot(ωi∇ωj−ωj∇ωi)=∇×(ωi∇ωj−ωj∇ωi)=∇ωi×∇ωj+ωi∇×∇ωj−(∇ωj×∇ωi+ωj∇×∇ωi)=∇ωi×∇ωj+ωirot(grad(ωj))−(∇ωj×∇ωi+ωjrot(grad(ωi))=2∇ωi×∇ωj\begin{aligned} \mathrm{rot}\bm{N}_ k &= \mathrm{rot}(\omega_i \nabla \omega_j - \omega_j \nabla \omega_i) \\ &= \nabla \times (\omega_i \nabla \omega_j - \omega_j \nabla \omega_i) \\ &= \nabla \omega_i \times \nabla \omega_j + \omega_i \nabla \times \nabla \omega_j \\ &-(\nabla \omega_j \times \nabla \omega_i + \omega_j \nabla \times \nabla \omega_i) \\ &= \nabla \omega_i \times \nabla \omega_j + \omega_i \mathrm{rot}(\mathrm{grad}(\omega_j)) \\ &-(\nabla \omega_j \times \nabla \omega_i + \omega_j \mathrm{rot}(\mathrm{grad}(\omega_i)) \\ &= 2 \nabla \omega_i \times \nabla \omega_j \end{aligned}

rotgrad=0\mathrm{rot}\mathrm{grad} = 0となることに気をつけると良い。 結果的に

rotNk=2∇ωi×∇ωj\mathrm{rot}\bm{N}_ k = 2 \nabla \omega_i \times \nabla \omega_j

座標の変換と全体行列

有限要素法では必ず座標の変換が必要になる。

(∂∂ξ1∂∂ξ2∂∂ξ3)=(∂ξ1∂ξ2∂ξ3),(∂∂x∂∂y∂∂z)=(∂x∂y∂z)\begin{pmatrix} \frac{\partial}{\partial\xi_1} \\ \frac{\partial}{\partial\xi_2} \\ \frac{\partial}{\partial\xi_3} \\ \end{pmatrix} = \begin{pmatrix} \partial_{\xi_1} \\ \partial_{\xi_2} \\ \partial_{\xi_3} \\ \end{pmatrix} , \begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \\ \end{pmatrix} = \begin{pmatrix} \partial_{x} \\ \partial_{y} \\ \partial_{z} \\ \end{pmatrix}

とする。

(∂ξ1∂ξ2∂ξ3)=(∂x∂ξ1∂y∂ξ1∂z∂ξ1∂x∂ξ2∂y∂ξ2∂z∂ξ2∂x∂ξ3∂y∂ξ3∂z∂ξ3)(∂x∂y∂z)\begin{pmatrix} \partial_{\xi_1} \\ \partial_{\xi_2} \\ \partial_{\xi_3} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial x}{\partial \xi_1} & \frac{\partial y}{\partial \xi_1} & \frac{\partial z}{\partial \xi_1} \\ \frac{\partial x}{\partial \xi_2} & \frac{\partial y}{\partial \xi_2} & \frac{\partial z}{\partial \xi_2} \\ \frac{\partial x}{\partial \xi_3} & \frac{\partial y}{\partial \xi_3} & \frac{\partial z}{\partial \xi_3} \\ \end{pmatrix} \begin{pmatrix} \partial_{x} \\ \partial_{y} \\ \partial_{z} \\ \end{pmatrix}

ここで、行列

J=(∂x∂ξ1∂y∂ξ1∂z∂ξ1∂x∂ξ2∂y∂ξ2∂z∂ξ2∂x∂ξ3∂y∂ξ3∂z∂ξ3)J = \begin{pmatrix} \frac{\partial x}{\partial \xi_1} & \frac{\partial y}{\partial \xi_1} & \frac{\partial z}{\partial \xi_1} \\ \frac{\partial x}{\partial \xi_2} & \frac{\partial y}{\partial \xi_2} & \frac{\partial z}{\partial \xi_2} \\ \frac{\partial x}{\partial \xi_3} & \frac{\partial y}{\partial \xi_3} & \frac{\partial z}{\partial \xi_3} \\ \end{pmatrix}

とする。

全体行列

有限要素法では以下のような行列を計算することがある。

rotNe⋅rotNe′Ne⋅Ne′∇ωi⋅Ne∇ωi⋅∇ωj\begin{aligned} & \mathrm{rot}\bm{N}_ e \cdot \mathrm{rot}\bm{N}_ {e'} \\ & \bm{N}_ e \cdot \bm{N}_ {e'} \\ & \nabla \omega_i \cdot \bm{N}_ e \\ & \nabla \omega_i \cdot \nabla \omega_j \\ \end{aligned}

ただし、以上でのrot,∇\mathrm{rot}, \nablaは全て(x,y,z)(x, y, z)においての微分である。 これを予め計算しやすくしておこう。 そのためには座標変換に注意する必要がある。 先に結果を示すと、

rotNe⋅rotNe′=4{⟨∇ξωi⋅∇ξωi′⟩⟨∇ξωj⋅∇ξωj′⟩−⟨∇ξωi⋅∇ξωj′⟩⟨∇ξωj⋅∇ξωi′⟩}Nex⋅Ne′x=⟨Neξ⋅Ne′ξ⟩∇ωi⋅Ne=⟨∇ξωi⋅Neξ⟩∇ωi⋅∇ωj=⟨∇ξωi⋅∇ξωj⟩\begin{align} \mathrm{rot}\bm{N}_ e \cdot \mathrm{rot}\bm{N}_ {e'} &= 4\{\langle\nabla_{\xi} \omega_i \cdot \nabla_{\xi} \omega_{i'}\rangle\langle\nabla_{\xi} \omega_j \cdot \nabla_{\xi} \omega_{j'}\rangle \nonumber \\ &- \langle\nabla_{\xi} \omega_i \cdot \nabla_{\xi} \omega_{j'}\rangle\langle\nabla_{\xi} \omega_j \cdot \nabla_{\xi} \omega_{i'}\rangle \} \\ \bm{N}^{x}_ e \cdot \bm{N}^{x}_ {e'} &= \langle \bm{N}^{\xi}_ e \cdot \bm{N}^{\xi}_ {e'} \rangle \\ \nabla \omega_i \cdot \bm{N}_ e &= \langle\nabla_{\xi}\omega_i\cdot \bm{N}^{\xi}_ e\rangle \\ \nabla \omega_i \cdot \nabla \omega_j &= \langle\nabla_{\xi}\omega_i\cdot\nabla_{\xi}\omega_j\rangle \\ \end{align}

となる。

rotNe⋅rotNe′\mathrm{rot}\bm{N}_ e \cdot \mathrm{rot}\bm{N}_ {e'}

rotNe=2∇ωi×∇ωj\mathrm{rot}\bm{N}_ {e} = 2\nabla \omega_i \times \nabla \omega_j

のように書けることは既に述べた。 ここで、ベクトル4重積を紹介する。

(A×B)⋅(C×D)=(A⋅C)(B⋅D)−(A⋅D)(B⋅C)(A\times B)\cdot(C\times D) = (A\cdot C)(B\cdot D) - (A\cdot D)(B\cdot C)

これを用いて、rotNe⋅rotNe′\mathrm{rot}\bm{N}_ e \cdot \mathrm{rot}\bm{N}_ {e'}を計算する。

rotNe⋅rotNe′=2(∇ωi×∇ωj)⋅2(∇ωi′×∇ωj′)=4(∇ωi×∇ωj)⋅(∇ωi′×∇ωj′)=4{(∇ωi⋅∇ωi′)(∇ωj⋅∇ωj′)−(∇ωi⋅∇ωj′)(∇ωj⋅∇ωi′)}=4{(∇ωi⋅∇ωi′)(∇ωj⋅∇ωj′)−(∇ωi⋅∇ωj′)(∇ωj⋅∇ωi′)}\begin{aligned} \mathrm{rot}\bm{N}_ e \cdot \mathrm{rot}\bm{N}_ {e'} &= 2(\nabla \omega_i \times \nabla \omega_j)\cdot 2(\nabla \omega_{i'} \times \nabla \omega_{j'}) \\ &= 4(\nabla \omega_i \times \nabla \omega_j)\cdot (\nabla \omega_{i'} \times \nabla \omega_{j'}) \\ &= 4\{(\nabla \omega_i \cdot \nabla \omega_{i'})(\nabla \omega_j \cdot \nabla \omega_{j'}) - (\nabla \omega_i \cdot \nabla \omega_{j'})(\nabla \omega_j \cdot \nabla \omega_{i'}) \} \\ &= 4\{(\nabla \omega_i \cdot \nabla \omega_{i'})(\nabla \omega_j \cdot \nabla \omega_{j'}) - (\nabla \omega_i \cdot \nabla \omega_{j'})(\nabla \omega_j \cdot \nabla \omega_{i'}) \} \\ \end{aligned}

ここで、係数を無視した以下を考える。また、∇=J−1∇ξ\nabla = J^{-1}\nabla_{\xi}であり、

∇ωi⋅∇ωj=(∇ωi)T(∇ωj)=(J−1∇ξωi)T(J−1∇ξωj)=(∇ξωi)T(J−TJ−1)∇ξωj=(∇ξωi)T(JJT)−1∇ξωj=(∇ξωi)Tg−1∇ξωj\begin{aligned} \nabla\omega_i \cdot \nabla \omega_j &= (\nabla\omega_i)^T(\nabla \omega_j) \\ & = (J^{-1}\nabla_{\xi}\omega_i)^T(J^{-1}\nabla_{\xi}\omega_j) \\ & = (\nabla_{\xi}\omega_i)^T(J^{-T}J^{-1})\nabla_{\xi}\omega_j \\ & = (\nabla_{\xi}\omega_i)^T(JJ^{T})^{-1}\nabla_{\xi}\omega_j \\ & = (\nabla_{\xi}\omega_i)^Tg^{-1}\nabla_{\xi}\omega_j \\ \end{aligned}

ここで、g=JJTg = JJ^{T}とする。 である。 また、⟨a⋅b⟩=aTg−1b\langle \bm{a} \cdot \bm{b} \rangle = \bm{a}^{T}g^{-1}\bm{b}とする。

(∇ωi×∇ωj)⋅(∇ωi′×∇ωj′)=(∇ωi⋅∇ωi′)(∇ωj⋅∇ωj′)−(∇ωi⋅∇ωj′)(∇ωj⋅∇ωi′)=⟨∇ξωi⋅∇ξωi′⟩⟨∇ξωj⋅∇ξωj′⟩−⟨∇ξωi⋅∇ξωj′⟩⟨∇ξωj⋅∇ξωi′⟩\begin{aligned} (\nabla \omega_i \times \nabla \omega_j)\cdot (\nabla \omega_{i'} \times \nabla \omega_{j'}) &= (\nabla \omega_i \cdot \nabla \omega_{i'})(\nabla \omega_j \cdot \nabla \omega_{j'}) \\ &- (\nabla \omega_i \cdot \nabla \omega_{j'})(\nabla \omega_j \cdot \nabla \omega_{i'}) \\ &= \langle\nabla_{\xi} \omega_i \cdot \nabla_{\xi} \omega_{i'}\rangle\langle\nabla_{\xi} \omega_j \cdot \nabla_{\xi} \omega_{j'}\rangle \\ &- \langle\nabla_{\xi} \omega_i \cdot \nabla_{\xi} \omega_{j'}\rangle\langle\nabla_{\xi} \omega_j \cdot \nabla_{\xi} \omega_{i'}\rangle \\ \end{aligned}

よって、

rotNe⋅rotNe′=4{⟨∇ξωi⋅∇ξωi′⟩⟨∇ξωj⋅∇ξωj′⟩−⟨∇ξωi⋅∇ξωj′⟩⟨∇ξωj⋅∇ξωi′⟩}\begin{aligned} \mathrm{rot}\bm{N}_ e \cdot \mathrm{rot}\bm{N}_ {e'} &= 4\{\langle\nabla_{\xi} \omega_i \cdot \nabla_{\xi} \omega_{i'}\rangle\langle\nabla_{\xi} \omega_j \cdot \nabla_{\xi} \omega_{j'}\rangle \\ &- \langle\nabla_{\xi} \omega_i \cdot \nabla_{\xi} \omega_{j'}\rangle\langle\nabla_{\xi} \omega_j \cdot \nabla_{\xi} \omega_{i'}\rangle \}\\ \end{aligned}

となる。

Ne⋅Ne′\bm{N}_ e \cdot \bm{N}_ {e'}

Nx=ωi∇ωj−ωj∇ωiNξ=ωi∇ξωj−ωj∇ξωiNx=J−1Nξ\begin{aligned} \bm{N}^{x} &= \omega_i \nabla \omega_j - \omega_j \nabla \omega_i \\ \bm{N}^{\xi} &= \omega_i \nabla_{\xi} \omega_j - \omega_j \nabla_{\xi} \omega_i \\ \bm{N}^{x} &= J^{-1}\bm{N}^{\xi} \end{aligned}

であるとする。

Nex⋅Ne′x=(Nex)TNe′x=(J−1Neξ)TJ−1Ne′ξ=(Neξ)TJ−TJ−1Ne′ξ=(Neξ)T(JJT)−1Ne′ξ=(Neξ)Tg−1Ne′ξ=⟨Neξ⋅Ne′ξ⟩\begin{aligned} \bm{N}^{x}_ e \cdot \bm{N}^{x}_ {e'} &= (\bm{N}^{x}_ e)^{T}\bm{N}^{x}_ {e'} \\ &= (J^{-1}\bm{N}^{\xi}_ e)^{T}J^{-1}\bm{N}^{\xi}_ {e'} \\ &= (\bm{N}^{\xi}_ e)^{T}J^{-T}J^{-1}\bm{N}^{\xi}_ {e'} \\ &= (\bm{N}^{\xi}_ e)^{T}(JJ^{T})^{-1}\bm{N}^{\xi}_ {e'} \\ &= (\bm{N}^{\xi}_ e)^{T}g^{-1}\bm{N}^{\xi}_ {e'} \\ &= \langle \bm{N}^{\xi}_ e \cdot \bm{N}^{\xi}_ {e'} \rangle \end{aligned}

となる。

∇ωi⋅Ne\nabla \omega_i \cdot \bm{N}_ e

∇ωi⋅Ne=J−1∇ξωi⋅J−1Neξ=(J−1∇ξωi)TJ−1Neξ=(∇ξωi)TJ−TJ−1Neξ=(∇ξωi)T(JJT)−1Neξ=(∇ξωi)Tg−1Neξ=⟨∇ξωi⋅Neξ⟩\begin{aligned} \nabla \omega_i \cdot \bm{N}_ e &= J^{-1}\nabla_{\xi}\omega_i \cdot J^{-1}\bm{N}^{\xi}_ e \\ &= (J^{-1}\nabla_{\xi}\omega_i)^{T}J^{-1}\bm{N}^{\xi}_ e \\ &= (\nabla_{\xi}\omega_i)^{T}J^{-T}J^{-1}\bm{N}^{\xi}_ e \\ &= (\nabla_{\xi}\omega_i)^{T}(JJ^{T})^{-1}\bm{N}^{\xi}_ e \\ &= (\nabla_{\xi}\omega_i)^{T}g^{-1}\bm{N}^{\xi}_ e \\ &= \langle\nabla_{\xi}\omega_i\cdot \bm{N}^{\xi}_ e\rangle \\ \end{aligned}

∇ωi⋅∇ωj\nabla \omega_i \cdot \nabla \omega_j

∇ωi⋅∇ωj=J−1∇ξωi⋅J−1∇ξωj=(J−1∇ξωi)TJ−1∇ξωj=(∇ξωi)TJ−TJ−1∇ξωj=(∇ξωi)T(JJT)−1∇ξωj=(∇ξωi)Tg−1∇ξωj=⟨∇ξωi⋅∇ξωj⟩\begin{aligned} \nabla \omega_i \cdot \nabla \omega_j &= J^{-1}\nabla_{\xi}\omega_i \cdot J^{-1}\nabla_{\xi}\omega_j \\ &= (J^{-1}\nabla_{\xi}\omega_i)^{T}J^{-1}\nabla_{\xi}\omega_j \\ &= (\nabla_{\xi}\omega_i)^{T}J^{-T}J^{-1}\nabla_{\xi}\omega_j \\ &= (\nabla_{\xi}\omega_i)^{T}(JJ^{T})^{-1}\nabla_{\xi}\omega_j \\ &= (\nabla_{\xi}\omega_i)^{T}g^{-1}\nabla_{\xi}\omega_j \\ &= \langle\nabla_{\xi}\omega_i\cdot\nabla_{\xi}\omega_j\rangle \\ \end{aligned}

有限要素行列

式(eq:A−ϕ)(eq: A-\phi)を行列にする。 物理法則はあらゆる空間のすべての点で成立する法則であるのでこの法則は要素mm内の任意の点x\bm{x}でも成り立つ。

rotμ−1rotA+σ∂A∂t+σgrad∂ϕ∂t=Jsource\mathrm{rot}\mu^{-1}\mathrm{rot}\bm{A} + \sigma\frac{\partial \bm{A}}{\partial t} + \sigma\mathrm{grad}\frac{\partial \phi}{\partial t} = \bm{J}_ {source} A(x)=∑e∈EdgeAeNe(x)\bm{A}(\bm{x}) = \sum_{e\in \mathrm{Edge}}A_e\bm{N}_ e(\bm{x}) ϕ(x)=∑n∈NodeϕiLi(x)\phi(\bm{x}) = \sum_{n \in {Node}}\phi_i L_i(\bm{x}) Nk=ωi∇ωj−ωj∇ωi\bm{N}_ k = \omega_i \nabla \omega_j - \omega_j \nabla \omega_i

この式にガラーキン法を適用する。 重みは辺Ne′\bm{N}_ {e'}だ。
ベクトル公式を一つ紹介する。 なお、時間による偏微分をみやすくするために∂A∂t=A˙\frac{\partial \bm{A}}{\partial t} = \dot{\bm{A}}とする。

div(a×b)=b⋅rota−a⋅rotb\mathrm{div}(\bm{a}\times \bm{b}) = \bm{b}\cdot\mathrm{rot}\bm{a} - \bm{a}\cdot\mathrm{rot}\bm{b}

それじゃlet’s go とりあえず両辺にNe′\bm{N}_ {e'}との内積を取る。

Ne′⋅(rotμ−1rotA+σA˙+σgradϕ˙)=Ne′⋅Jsource\bm{N}_ {e'}\cdot (\mathrm{rot}\mu^{-1}\mathrm{rot}\bm{A} + \sigma\dot{\bm{A}} + \sigma\mathrm{grad}\dot{\phi} )= \bm{N}_ {e'}\cdot\bm{J}_ {source}

右辺から計算していこう

(LHS)=μ−1rotA⋅rotNe′+div(μ−1rotA×Ne′)+σ(Ne′⋅A˙+gradϕ˙)\begin{aligned} (LHS) &= \mu^{-1}\mathrm{rot}\bm{A}\cdot\mathrm{rot}\bm{N}_ {e'} + \mathrm{div}(\mu^{-1}\mathrm{rot}\bm{A}\times\bm{N}_ {e'}) \\ &+\sigma(\bm{N}_ {e'}\cdot\dot{\bm{A}} + \mathrm{grad}\dot{\phi}) \\ \end{aligned}

これを要素mm内のx\bm{x}について積分(要素mm内、体積分)する。 これを解析対象MMについて積分(体積分)する。解析対象MMが占める三次元空間をΩ\Omegaとする。

∫ΩNe′⋅(rotμ−1rotA+σA˙+σgradϕ˙)dV=∫ΩNe′⋅JsourcedV\int_{\Omega}\bm{N}_ {e'}\cdot (\mathrm{rot}\mu^{-1}\mathrm{rot}\bm{A} + \sigma\dot{\bm{A}} + \sigma\mathrm{grad}\dot{\phi})\mathrm{d}V = \int_{\Omega}\bm{N}_ {e'}\cdot\bm{J}_ {source}\mathrm{d}V

先程計算したものを使用して右辺を計算する。

(LHS)=∫Ωμ−1rotA⋅rotNe′+div(μ−1rotA×Ne′)+σ(Ne′⋅A˙+Ne′⋅gradϕ˙)dV=∫Ωμ−1rotA⋅rotNe′dV+∫Ωdiv(μ−1rotA×Ne′)dV+∫Ωσ(Ne′⋅A˙+Ne′⋅gradϕ˙)dV=∫Ωμ−1rotA⋅rotNe′dV+∫∂Ωμ−1rotA×Ne′⋅dS+∫Ωσ(Ne′⋅A˙+Ne′⋅gradϕ˙)dV\begin{aligned} (LHS) &= \int_{\Omega}\mu^{-1}\mathrm{rot}\bm{A}\cdot\mathrm{rot}\bm{N}_ {e'} + \mathrm{div}(\mu^{-1}\mathrm{rot}\bm{A}\times\bm{N}_ {e'}) \\ &+\sigma(\bm{N}_ {e'}\cdot\dot{\bm{A}} + \bm{N}_ {e'}\cdot\mathrm{grad}\dot{\phi}) \mathrm{d}V\\ &= \int_{\Omega}\mu^{-1}\mathrm{rot}\bm{A}\cdot\mathrm{rot}\bm{N}_ {e'}\mathrm{d}V + \int_{\Omega}\mathrm{div}(\mu^{-1}\mathrm{rot}\bm{A}\times\bm{N}_ {e'})\mathrm{d}V \\ &+ \int_{\Omega}\sigma(\bm{N}_ {e'}\cdot\dot{\bm{A}} + \bm{N}_ {e'}\cdot\mathrm{grad}\dot{\phi}) \mathrm{d}V \\ &= \int_{\Omega}\mu^{-1}\mathrm{rot}\bm{A}\cdot\mathrm{rot}\bm{N}_ {e'}\mathrm{d}V + \int_{\partial \Omega}\mu^{-1}\mathrm{rot}\bm{A}\times\bm{N}_ {e'}\cdot\mathrm{d}\bm{S} \\ &+ \int_{\Omega}\sigma(\bm{N}_ {e'}\cdot\dot{\bm{A}} + \bm{N}_ {e'}\cdot\mathrm{grad}\dot{\phi}) \mathrm{d}V \\ \end{aligned}

ここで二段目の式変形の第二項ではガウスの法則を使った。 境界では磁束密度B=rotA\bm{B}=\mathrm{rot}\bm{A}は面に対して垂直であり特別に∂Ω\partial\OmegaにおいてB∥n\bm{B}\parallel\bm{n}つまりB×n=0\bm{B}\times \bm{n}=0 という条件をつける。 これにより、

(LHS)=∫Ωμ−1rotA⋅rotNe′dV+∫Ωσ(Ne′⋅A˙+Ne′⋅gradϕ˙)dV\begin{aligned} (LHS) &= \int_{\Omega}\mu^{-1}\mathrm{rot}\bm{A}\cdot\mathrm{rot}\bm{N}_ {e'}\mathrm{d}V + \int_{\Omega}\sigma(\bm{N}_ {e'}\cdot\dot{\bm{A}} + \bm{N}_ {e'}\cdot\mathrm{grad}\dot{\phi}) \mathrm{d}V \\ \end{aligned}

さて、積分の前提として積分対象は非交叉な積分対象に分割できる。

Ω=∑mVm\Omega = \sum_{m}V_m

とする。
ここで分割した積分対象 要素mmについて考える。 ここでA,ϕ\bm{A}, \phiなどを形状関数で表現したものを代入する。なお、第二項と第三項を入れ替えている。

(LHS)m=∫Vmμ−1rotA⋅rotNe′dV+∫Vmσ(Ne′⋅A˙+Ne′⋅gradϕ˙)dV=∑e∈mEAe∫Vmμ−1rotNe⋅rotNe′dV+∑e∈mEA˙e∫VmσNe′⋅NedV+∑n∈mNϕ˙n∫VmσNe′⋅gradLndV\begin{aligned} (LHS)_m &= \int_{V_m}\mu^{-1}\mathrm{rot}\bm{A}\cdot\mathrm{rot}\bm{N}_ {e'}\mathrm{d}V + \int_{V_m}\sigma(\bm{N}_ {e'}\cdot\dot{\bm{A}} + \bm{N}_ {e'}\cdot\mathrm{grad}\dot{\phi}) \mathrm{d}V \\ &= \sum_{e\in m_E} A_ e\int_{V_m}\mu^{-1}\mathrm{rot}\bm{N}_ e\cdot\mathrm{rot}\bm{N}_ {e'}\mathrm{d}V \\ &+ \sum_{e\in m_E}\dot{A}_ e\int_{V_m}\sigma\bm{N}_ {e'}\cdot\bm{N}_ e \mathrm{d}V + \sum_{n\in m_N}\dot{\phi}_ n\int_{V_m}\sigma\bm{N}_ {e'}\cdot\mathrm{grad}L_n \mathrm{d}V \\ \end{aligned}

そこで

Kee′=∫Vmμ−1rotNe⋅rotNe′dVKnn=0Ke′n=Kne′=0Ce′e=Cee′=∫VmσNe′⋅NedVCe′n=∫VmσNe′⋅gradLndV\begin{aligned} K_{ee'} &= \int_{V_m}\mu^{-1}\mathrm{rot}\bm{N}_ e\cdot\mathrm{rot}\bm{N}_ {e'}\mathrm{d}V \\ K_{nn} &= 0 \\ K_{e'n} = K_{ne'} &= 0 \\ C_{e'e} = C_{ee'} &= \int_{V_m}\sigma\bm{N}_ {e'}\cdot\bm{N}_ {e}\mathrm{d}V \\ C_{e'n} &= \int_{V_m}\sigma\bm{N}_ {e'}\cdot\mathrm{grad}L_n\mathrm{d}V \end{aligned}

これを用いると

(LHS)=∑m(LHS)m(LHS)m=∑eAeKee′+∑eA˙eCe′e+∑nϕ˙nCe′n\begin{aligned} (LHS) &= \sum_m(LHS)_m \\ (LHS)_m &= \sum_e A_eK_{ee'} + \sum_e\dot{A}_ eC_{e'e} + \sum_{n}\dot{\phi}_ nC_{e'n} \end{aligned}

また、Kee′,Ke′n,Knn,Cee′,Ce′n,CnnK_{ee'},K_{e'n},K_{nn},C_{ee'},C_{e'n},C_{nn}などを成分にもつ行列を考え、

A=(A0,…,Ae,…,AE−1,ϕ0,…,ϕn…ϕN−1)\mathbf{A}=(A_0,\ldots, A_e,\ldots, A_{E-1}, \phi_0, \ldots, \phi_n\ldots\phi_{N-1})

とし、

(RHS)=∫ΩNe′⋅JsourcedV=Fe′\begin{aligned} (RHS) &= \int_{\Omega}\bm{N}_ {e'}\cdot\bm{J}_ {source}\mathrm{d}V \\ &= F_e' \end{aligned}

とする。

KA+CA˙=F\mathrm{K}\mathbf{A} + \mathrm{C}\dot{\mathbf{A}} = \mathrm{F}

が得られる。

不完全コレスキー分解前処理共役勾配法

ゲージ変換と電流の発散条件

さて、先程導出した式だが実はうまく解けない。 時間変化をしないC=O\mathrm{C}=\mathrm{O}の状態でも解けない。 これはK\mathrm{K}が必ず特異行列になることに由来する。 これは四面体を考えてみればわかる。

対象から四面体要素mmを一つ取り出して磁束密度B=rotA\bm{B}=\mathrm{rot}\bm{A}を計算することを考えてみよう。 辺形状関数の標本点は辺に割当てられ、磁束密度は面に割り当てれられる。 磁束密度がベクトルポテンシャルの標本点の線形和で表されるとする。 四面体の面は44つに対して辺は66本だ。つまり、各面に割り当てられる磁束密度を求めるにあたって辺の数が多すぎるため方程式は不定となる。(自由度が多すぎて解が一意に定まらない。)

四面体のある面でループを作る際、それは補木と木に分割することができる。木(または補木)を一意に決定するとループが決まるため、KKのrank\mathrm{rank}は補木の辺の数e−n+1e-n+1と等しくなる。 これは図形の根本的な問題をはらんでいる。

そこで以下ではどのようにすればこの方程式が解けるのかを見ていく。

クランク・ニコルソン法

共役勾配法

不完全コレスキー分解