超伝導体のさまざまな電磁現象のFEM による解析

木内研究室

内倉 信聖

平成 25 年 2 月 25 日

電子情報工学科

目次
1

第 1 章 序論
1.1

超伝導体 ........................................................................................................................ 1

1.2

有限要素法(FEM : Finite Element Method).............................................................. 1

1.3

超伝導体内の電磁現象の数値解析の方法..................................................................... 2

1.4

メッシュサイズの注意点 .............................................................................................. 41.5

弱形式 ........................................................................................................................... 4

1.6

支配方程式の弱形式化 .................................................................................................. 5

1.7

ガラーキン法 ................................................................................................................ 81.7.1

近似解 .................................................................................................................... 81.7.2

重み付き残差法 ..................................................................................................... 9

1.7.3

ガラーキン法 ......................................................................................................... 9

1.8

グリーンの法則 ............................................................................................................. 9

1.8.1

ガウスの発散定理 .................................................................................................. 9

1.8.2

グリーンの法則の証明......................................................................................... 101.9磁束密度 B の計算 .......................................................................................................11

1.10

PHOTO-Series ........................................................................................................... 12

1.11

本研究の目的 .............................................................................................................. 1213

第 2 章 解析方法
2.1

円筒形超伝導体に電流を軸方向に流す場合 ............................................................... 13

2.2

コイルに電流を流す場合 ............................................................................................ 1619

第 3 章 結果と考察
3.1

2.1 のモデルにおける磁束密度 B の分布 ................................................................... 19

3.2

2.2 のモデルにおける磁束密度 B の分布 ................................................................... 2326

第 4 章 まとめ

i

表目次

ii

図目次
1.1 有限要素法の概念 ........................................................................................................... 2
1.2 高温超伝導体における I-V 特性と数値解析コードでの低電流領域の近似特性 ............. 3
2.1 円筒形超伝導体の全体図と解析モデル......................................................................... 132.2 図 2.1 の解析モデルの寸法 ........................................................................................... 14
2.3 PHOTO-Series を用いて作製した円筒超伝導体の解析モデル.................................... 14
2.4 電流𝐼を流した時の磁束密度𝐵の分布の様子 ................................................................. 15
2.5 超伝導体に流れる交流電流のsin𝜔𝑡の波形 ................................................................... 15
2.6 超伝導体の上にコイルがある状態のモデル ................................................................. 162.7 図 2.6 の解析モデルの寸法 ........................................................................................... 17
2.8 PHOTO-Series を用いて作製した解析モデル ............................................................. 17
2.9 時間毎のコイルと超伝導体の距離 ................................................................................ 18
3.1 PHOTO-Series を用いて作製したモデルで円筒超伝導体の磁束密度ベクトルの様子 19
3.2 円筒超伝導体内の時間毎の磁束密度 ............................................................................ 20
3.3 𝑛 = 20, 200の時の磁束密度 ........................................................................................... 21
3.4 𝜎max⁄𝜎init = 10, 100, 1000の時の磁束密度 .................................................................... 22
3.5 PHOTO-Series を用いて作製したモデルで超伝導体にコイルを近づけた時の磁束密度ベクトルの様子............................................................................................... 23
3.6 コイルを超伝導体に近づけなかった時の超伝導体内の磁束密度................................. 24
3.7 コイルを超伝導体に近づけた時の超伝導体内の磁束密度 ........................................... 25

iii

第1章
1.1

序論

超伝導体

1908 年、オランダのヘイケ・カメルリング・オネス (Heike Kamerlingh Onnes) がヘリ
ウムの液化に初めて成功した．さらに、液体ヘリウムが極低温であることにより、1911 年に水銀の抵抗が 4 K 付近で突然ゼロになる現象を発見した．このような物理現象は今までにないということが分かり、この現象を持つ物質は超伝導体と呼ばれるようになった．超伝導状態では電気抵抗が無いことで大電流が通電できるということが期待され、しばらくの間研究が行われたが、決定的な理論は発見されなかった．しかし、BCS 理論が 1957 年に発見されたことにより、超伝導体の発見機構が明らかにされた．この理論では超伝導体が超伝導状態から常伝導状態へと転移する温度である臨界温度𝑇c は 30 K 程度だと考えられた．しかしながら、1986 年に 30 K 以上の𝑇c を持つ La2-xBaxCuO4 などの La-Ba-Cu-O 系超伝導体が、ドイツの物理学者・鉱物学者のヨハネス・ゲオルク・ベドノルツ (Johannees Georg Bednorz) とスイスの物理学者のカール・アレクサンダー・ミュラー (Karl Alexander Mü ller) によって発見された．その後𝑇c が液体窒素の沸点(77.3 K)以上の高温超伝導体も発見され、今後より高い𝑇c を持つ超伝導体の発見が期待されている．

1.2

有限要素法(FEM : Finite Element Method)

有限要素法(FEM：Finite Element Method)は、解析的に解くことが難しい微分方程式の近似解を数値的に得る方法の１つである．円柱や無限平板のような単純な形状ではなく、複雑な形状の問題だと解析的に解くことは非常に困難である．そこで複雑な形状の問題の解析を行う場合は、対象物を単純な形状の要素の集合体としてとらえ、各々の要素間で境界条件を満たすように方程式を作製する．それから、それぞれの要素で作製された方程式を対象物全体の連
立一次方程式として組み立てて計算を行う．またメッシュと呼ばれる分割された要素を細かくすることで、計算精度は増加する．しかし、有限要素法は単なる数値解析手法であるため、解析対象物のモデリングが適切でないと間違った解析結果を導く可能性が高い．そのため、解析対象物についてよく理解しておく必要がある．有限要素法の概念を図 1.1 に示す．

1

複雑な図形

分割それぞれ計算再構成

図 1.1: 有限要素法の概念

1.3

超伝導体内の電磁現象の数値解析の方法

超伝導体の抵抗は、電流密度が臨界電流密度に近づくにつれてゆるやかに発生する特性を持つ．このため、電流-電圧特性において電界が電流密度の𝑛乗に比例するとした𝑛値モデルを仮定する．電磁界に関する支配方程式は、透磁率を𝜇、導電率を𝜎とするとベクトルポテンシャ
ル A とスカラーポテンシャル𝜙を用いて次式で表すことができる．
1∇ × ∇ × 𝑨 = −𝜎(𝑨̇ + ∇𝜙)
𝜇

(1.1)

∇ ⋅ 𝜎(𝑨̇ + ∇𝜙) = 0(1.2)

また、磁界 B、電界 E、電流密度 J はそれぞれ次式から導出される．
𝑩 = ∇ × 𝑨、𝑬 = −𝛻𝜙 − 𝑨̇、𝑱 = 𝜎𝑬

(1.3)

電界基準𝐸0 を発生している時の電流密度を𝐽cとすると電流-電圧特性は、以下のように表せる．
𝐽 𝑛
𝐸 = 𝐸0 ( )
𝐽c

(1.4)

指数𝑛は超伝導から常伝導に転移するときの鋭さを表す．𝑛値が小さい場合、電流値を少し下げた時に発生する電圧は急速に小さくならない．しかし𝑛値が大きい場合は、電流値を少し下げた時に発生する電圧は急速に小さくなる．このため、𝑛値が大きいほうが一般的に超伝導体の応用にとって良いとされている．
ここで、この強い非線形性を扱うために式(1.1)、式(1.2)、式(1.3)の𝜎の代わりに仮想的な
電気伝導率
𝜎s =

𝐽
𝐽c 𝐽c 𝑛−1
= ( )𝐸 𝐸0 𝐽

(1.5)

𝐽c
𝐸0

(1.6)

を用いる．ここで初期電気伝導率は
𝜎init =とし、また電気伝導率の上限値として最大電気伝導率を𝜎maxと定義する．

2

𝜎s は非線形性を持つので、上記に基づいて各解析時刻において収束反復演算を行う．また、その反復演算手順を以下に示す．
1. 初期伝導率𝜎init を全要素に設定
2. 電磁界方程式を解く
3. σs を計算し、前回の𝜎s との差を計算する
4. 差が小さくなるまで収束反復演算
5. 𝜎s が𝜎maxを超える場合は𝜎maxに置換する
収束の条件を以下に示す．
∑∆𝜎
< 𝜀
∑𝜎

(1.7)

ただし、𝜀はあらかじめ設定した基準値である．
ここで、初期伝導率について考えると、式(1.6)において𝐸0 または𝜎init を変えれば𝐽cの値も変化する．この現象はブロードな I-V 特性を持つ高温超伝導体に特に顕著に表れる．高温超伝導体の I-V 特性は以下の式で表される．
𝐽 𝑛𝐽c 𝐽 𝑛𝐸 = 𝐸0 ( ) =( )𝐽c𝜎init 𝐽c

(1.8)

また、超伝導体の I-V (E-J )特性と解析時の近似ラインを図 1.2 に示す．図 1.2 から低電流域において、近似ラインと I-V 特性カーブにずれがあり、結果的に低電流域での交流損失を多めに計算することになる．しかし、𝜎maxを大きく設定することでこの問題を解決できる．

図 1.2: 高温超伝導体における I-V 特性と数値解析コードでの低電流領域の近似特性(本庄卓一、「超電導ケーブルに生じる交流損失の数値解析」(JMAG Users Conference 2001)、
p.6‐5)
3

1.4

メッシュサイズの注意点

導体内での急激な電磁界分布の変動や、及びそれに伴う損失を精度良く評価するためには、モデルのメッシュサイズを細かく設定する必要がある．少なくとも以下の式で表される表皮厚さ𝛿以下の大きさに設定することが求められる．
2
2
𝛿=√
=√
𝜔𝜎𝜇
2𝜋𝑓𝜎𝜇

(1.9)

ただし、𝜎は電気伝導率、𝜇は透磁率である．
ここで表皮厚さとは、高周波電流が導体内を流れるときに、電流密度が導体表面で高く離れると低くなるという表皮効果において、導体の電流密度 J が表面からの深さ𝑥に対して以下のように減少する𝑥(1.10)
𝐽 = 𝑒 −𝐷における𝛿のことである．ここでの𝑥は電流が表面電流の 1/e(約 0.37)となる深さである．また eは自然数とする．
また超伝導体のモデルの要素分割において、上記の式中の𝜎を最大伝導率𝜎maxで置き換えた
値がメッシュサイズとなる．さらに超伝導体において𝜇は一般的に𝜇0 を用いる．

1.5

弱形式

有限要素法では微分方程式を後述するグリーンの法則を用いて、2 階の導関数を含む式を 1階までの導関数のみの形に変形する．これを弱形式化という．そのため、ここでは弱形式について説明する．弱形式とは微分方程式の導関数が低い方程式のことである．例として以下の条件を満たす𝑢を求める問題を示す．
強形式の場合−d2 𝑢= 𝑓(𝑥)
d𝑥 2(1.11)(0 < 𝑥 < 𝑎)(1.12)

𝑢(0) = 𝑢0 , 𝑢(𝑎) = 𝑢𝑎
弱形式の場合𝑎𝑎d𝑢 d𝑣
d𝑥 = ∫ 𝑓(𝑥)𝑣d𝑥
0 d𝑥 d𝑥0∫∀𝑣(1.13)(1.14)

𝑢(0) = 𝑢(𝑎) = 0
ただし、𝑣は𝑣(0) = 𝑣(𝑎) = 0を満たす 2 階微分可能な関数．

強形式を弱形式に変換する場合、強形式の両辺に𝑣(0) = 𝑣(𝑎) = 0を満たす 2 階微分可能な関数𝑣をかけ定義域全体(0から𝑎)まで部分積分する．
𝑎
d𝑢 d𝑣
d𝑢 𝑎
d𝑥 − [ 𝑣] = ∫ 𝑓(𝑥)𝑣d𝑥
d𝑥 0
0 d𝑥 d𝑥
0
𝑎

∫

4

∀𝑣

(1.15)

𝑣(0) = 𝑣(𝑎) = 0から左辺第 2 項は 0 なので、
𝑎

𝑎
d𝑢 d𝑣
d𝑥 = ∫ 𝑓(𝑥)𝑣d𝑥
0 d𝑥 d𝑥0∫∀𝑣(1.16)

上記の式は𝑣(0) = 𝑣(𝑎) = 0の境界条件さえ満たせば任意の𝑣について成立する．よって強形式
から弱形式に変換することができる．

1.6

支配方程式の弱形式化

有限要素法を用いて解析を行う為に必要な支配方程式の弱形式化する手順を以下に示す．
円柱座標系におけるスカラーポテンシャル𝜙の勾配(gradient)、およびベクトルポテンシャル Aの発散(divergence)、回転(rotation)はそれぞれ以下のように定義される．
𝜕𝜙
1 𝜕𝜙
𝜕𝜃
𝑒𝑟 +
𝑒𝜃 +
𝑒
𝜕𝑟
𝑟 𝜕𝜃
𝜕𝑧 𝑧
𝐴𝑟
1𝜕
1 𝜕𝐴𝜃 𝜕𝐴𝑧
1 𝜕 𝜕𝜙(𝑟𝐴𝑟 ) ++] [𝐴𝜃 ] =
𝑟 𝜕𝑟
𝑟 𝜕𝜃
𝜕𝑧
𝑟 𝜕𝜃 𝜕𝑧 𝐴
𝑧
∇𝜙 =

𝜕
∇⋅𝑨=[
𝜕𝑟

𝑒𝑟
𝜕
∇×𝑨= |
𝜕𝑟
𝐴𝑟

𝑒𝜃
1 𝜕
𝑟 𝜕𝜃
𝐴𝜃

1 𝜕𝐴𝑧 𝜕𝐴𝜃
{
−
}𝑒
𝑒𝑧
𝑟 𝜕𝜃
𝜕𝑧 𝑟
𝜕𝜙
𝜕𝐴𝑟 𝜕𝐴𝑧
|=
{
−
}𝑒
𝜕𝑧
𝜕𝑧
𝜕𝑟 𝜃
𝐴𝑧
1 𝜕
𝜕𝐴𝑟
{ (𝑟𝐴𝜃 ) −}𝑒[ 𝑟 𝜕𝑟𝜕𝜃 𝑧 ](1.17)(1.18)(1.19)

ここで𝑒r , 𝑒θ , 𝑒z はそれぞれ𝑟, 𝜃, 𝑧方向の単位ベクトルとする．また軸対称場ではベクトルポ
テンシャル A は𝜃方向成分のみを有することになり、物理量は𝜃方向には一様である．つまり以下のような仮定が成り立つ．
𝐴𝑟 = 𝐴𝑧 = 0、

∂𝐴𝜃
=0
𝜕𝜃

(1.20)

式(1.20)を式(1.19)に代入すると次のようになる．
𝜕𝐴𝜃
𝑒
𝜕𝑧 𝑟
∇×𝑨=
0
1𝜕
(𝑟𝐴𝜃 )𝑒𝑧[𝑟 𝜕𝑟](1.21)

ここで、式(1.20)を考慮して式(1.1)を計算すると、次式のように表される．
𝑒𝑟
𝑒𝜃
𝑒𝑧
𝜕
1 𝜕
𝜕𝜙
1
|
|
𝜕𝑟
𝑟 𝜕𝜃
𝜕𝑧
∇× ∇×𝑨=| 1 1 𝜕𝐴𝜇𝜕𝐴𝜃
1 𝜕𝐴𝑟 𝜕𝐴𝑧
1 1 𝜕
𝜕𝐴𝑟 |
𝑧
(
−
)
(
−
)
}]
[ { (𝑟𝐴𝜃 ) −𝜇 𝑟 𝜕𝜃𝜕𝑧
𝜇 𝜕𝑧
𝜕𝑟
𝜇 𝑟 𝜕𝑟𝜕𝜃
0
𝜕 1 𝜕𝐴𝜃𝜕 1 𝜕(𝑟𝐴𝜃 )}} 𝑒𝜃= {− ()− {𝜕𝑧 𝜇 𝜕𝑧
𝜕𝑟 𝜇𝑟 𝜕𝑟[]05

(1.22)

以上より、式(1.1)の具体的な成分を求めることができた．また、以降は𝐴𝜃 は特に断らない
限り A と表記する．式(1.1)に式(1.3)、(1.22)を代入して整理すると、支配方程式は次のように表される．
∂ 1 𝜕𝐴𝜕 1 𝜕(𝑟𝐴)} − 𝜎(𝐴̇ + ∇𝜙) = 0()+ {∂z 𝜇 𝜕𝑧
𝜕𝑟 𝜇𝑟 𝜕𝑟

(1.23)

ここで、要素の重み関数を{N }とすると、ガラーキン法を用いた式は以下のように表される．また、軸対称場であるので各項が2𝜋𝑟倍される．
∬{𝑁} [S𝜕 1 𝜕𝐴𝜕 1 𝜕(𝑟𝐴)}] 2𝜋𝑟d𝑟d𝑧()+ {𝜕𝑧 𝜇 𝜕𝑧
𝜕𝑟 𝜇𝑟 𝜕𝑟

(1.24)− ∬{𝑁}𝜎(𝐴̇ + ∇𝜙)2𝜋𝑟d𝑟d𝑧 = 0
S

ここで、式(1.24)の左辺第 1 項について考える．
𝜕 1 𝜕𝐴𝜕 1 𝜕(𝑟𝐴)}] 2𝜋𝑟d𝑟d𝑧∬{𝑁} [ (
)+ {
𝜕𝑧
𝜇
𝜕𝑧
𝜕𝑟
𝜇𝑟 𝜕𝑟
S

(1.25)

式(1.25)は 2 階微分の項があるので、弱形式化つまり 1 階微分形に次数を下げるために部
分積分を行う．部分積分の公式を導出するために次のスカラー関数 U と V の関係を示す公式
を考える．
∇ ⋅ (𝑈∇𝑉) = ∇𝑈 ⋅ ∇𝑉 + 𝑈 ⋅ ∇2 𝑉

(1.26)

この公式を次のように変形する．
∬ ∇ ⋅ (𝑈∇𝑉)d𝑆 = ∬(∇𝑈 ⋅ ∇𝑉 + 𝑈 ⋅ ∇2 𝑉)d𝑆
S

S

(1.27)

ここで式(1.27)の左辺においてグリーンの定理を考えると以下のようになる．
∬ ∇ ⋅ (𝑈∇𝑉)d𝑆 = ∮(𝑈∇𝑉) ⋅ 𝑛d𝑆
S(1.28)

C𝑛は有限要素領域の境界Γ上の法線ベクトルであり、その𝑟、𝑧方向成分を𝑛r 、𝑛z とする．式(1.28)を式(1.27)の左辺に適用し整理を行うと、以下の部分積分の公式が成り立つ．
∬(𝑈 ⋅ ∇2 𝑉)d𝑆 = ∮(𝑈∇𝑉) ⋅ 𝑛d𝑆 − ∬(∇𝑈 ⋅ ∇𝑉)d𝑆
S

CS(1.29)

上記の部分積分の公式おいて考える．
𝑈 = 2πr{𝑁}、𝜕 1 𝜕
𝜕 1 𝜕𝐴(𝑟𝐴)} + ({)
𝜕𝑟 𝜇𝑟 𝜕𝑟𝜕𝑧 𝜇 𝜕𝑧

(1.30)

1 𝐴 + 𝑟𝜕𝐴
[
𝜕𝑟 ]
𝜇𝑟
𝑟

(1.31)∇2 𝑉 =とすると、
𝑈∇𝑉 = (2𝜋𝑟{𝑁})(𝑈∇𝑉) ⋅ 𝑛 = (2𝜋𝑟{𝑁})

1
𝜕𝐴
[𝐴+𝑟
𝜇𝑟
𝜕𝑟

= (2πr{N}) [(

𝑟

𝑛𝑟
𝜕𝐴
] [𝑛 ]𝑧
𝜕𝑧

𝐴 1 𝜕𝐴𝜕𝐴+) 𝑛𝑟 +𝑛 ]𝜇𝑟 𝜇 𝜕𝑟𝜕𝑧 𝑧

6(1.32)

𝜕𝐴
1 𝐴 + 𝑟 𝜕𝑟𝜕𝜕
∇𝑈 ⋅ ∇𝑉 = [ (2𝜋𝑟{𝑁})](2𝜋𝑟{𝑁}) ] [
𝜕𝐴
𝜇𝑟
𝜕𝑟
𝜕𝑧
𝑟
𝜕𝑧
𝜕𝐴
𝜕{𝑁} {𝑁} 𝜕{𝑁} 1 𝐴 + 𝑟 𝜕𝑟
= 2𝜋𝑟 [
]
] [
+
𝜕𝐴
𝜕𝑟
𝑟
𝜕𝑧 𝜇𝑟
𝑟
𝜕𝑧
{𝑁}𝐴 {𝑁} 𝜕𝐴
𝜕{𝑁} 𝐴 1 𝜕𝐴
= 2𝜋𝑟 [
( +)+( 2 +)
𝜕𝑟 𝜇𝑟 𝜇 𝜕𝑟
𝜇𝑟
𝜇𝑟 𝜕𝑟
𝜕{𝑁} 1 𝜕𝐴+()]
𝜕𝑧 𝑟 𝜕𝑧

(1.33)となり、式(1.25)は以下のように表すことが出来る．
𝜕 1 𝜕𝐴𝜕 1 𝜕(𝑟𝐴)}] 2𝜋𝑟d𝑟d𝑧∬{𝑁} [ ()+ {𝜕𝑧 𝜇 𝜕𝑧
𝜕𝑟 𝜇𝑟 𝜕𝑟S
= ∮ 2𝜋𝑟{𝑁} [(C− ∬ 2𝜋𝑟 [S

𝐴 1 𝜕𝐴𝜕𝐴+) 𝑛𝑟 +𝑛 ] d𝑠𝜇𝑟 𝜇 𝜕𝑟𝜕𝑧 𝑧{𝑁}𝐴 {𝑁} 𝜕𝐴
𝜕{𝑁} 𝐴 1 𝜕𝐴
𝜕{𝑁} 1 𝜕𝐴( +)+( 2 +)+()] d𝑆
𝜕𝑟 𝜇𝑟 𝜇 𝜕𝑟
𝜇𝑟
𝜇𝑟 𝜕𝑟
𝜕𝑧 𝑟 𝜕𝑧

(1.34)

式(1.34)の右辺第 1 項は境界積分項である．ここで、自然境界条件は磁界 H が境界に対して垂直となる条件であり、以下のようになる．
𝑯×𝑛 =0(1.35)1 𝜕𝐴𝜇 𝜕𝑧
1𝑯= 𝛻×𝐴=
0
𝜇
1 𝜕
(𝑟𝐴)
[ 𝜇𝑟 𝜕𝑟
](1.36)

また、磁界 H は以下のように表せる．
−

ここで、H と𝒏の外積を考える．
𝑒𝑟
1 𝜕𝐴
𝑯×𝑛 = |−
𝜇 𝜕𝑧
𝑛𝑟
=[

𝑒𝜃
0
0

𝑒𝑧
1 𝜕
(𝑟𝐴) |
𝜇𝑟 𝜕𝑟
𝑛𝑧

1 𝜕
1 𝜕𝐴
(𝑟𝐴)𝑛𝑟 +𝑛 ]𝜇𝑟 𝜕𝑟𝜇 𝜕𝑧 𝑧

0
1
1 𝜕𝐴
1 𝜕𝐴
=[( 𝐴+) 𝑛𝑟 +𝑛 ]𝜇𝑟
𝜇 𝜕𝑟
𝜇 𝜕𝑧 𝑧0(1.37)

したがって、式(1.29)の条件より
𝐴 1 𝜕𝐴1 𝜕𝐴( +) 𝑛𝑟 +𝑛 =0
𝜇𝑟 𝜇 𝜕𝑟𝜇 𝜕𝑧 𝑧7(1.38)

これより、式(1.34)の右辺第 1 項の境界積分は0となる．有限要素法による磁場解析では、
この項は通常 0 となる．これによって以下の式が得られる．
𝜕 1 𝜕𝐴𝜕 1 𝜕(𝑟𝐴)}] 2𝜋𝑟d𝑟d𝑧∬{𝑁} [ (
)+ {
𝜕𝑧
𝜇
𝜕𝑧
𝜕𝑟
𝜇𝑟 𝜕𝑟
S
= − ∬ 2𝜋𝑟 [S{𝑁}𝐴 {𝑁} 𝜕𝐴
𝜕{𝑁} 𝐴 1 𝜕𝐴
𝜕{𝑁} 1 𝜕𝐴( +)+( 2 +)+()] d𝑆
𝜕𝑟 𝜇𝑟 𝜇 𝜕𝑟
𝜇𝑟
𝜇𝑟 𝜕𝑟
𝜕𝑧 𝑟 𝜕𝑧{𝑁}𝐴
𝐴 𝜕{𝑁}
𝑟 𝜕{𝑁} 𝜕𝐴d𝑟d𝑧 − 2𝜋 ∬d𝑟d𝑧 − 2𝜋 ∬
d𝑟d𝑧
𝜇
𝜕𝑟
𝜇
𝜕𝑟
𝜕𝑟
S
S
S 𝜇𝑟

= −2𝜋 ∬

{𝑁} 𝜕𝐴𝑟 𝜕{𝑁} 𝜕𝐴− 2𝜋 ∬
d𝑟d𝑧 − 2𝜋 ∬
d𝑟d𝑧
𝜇
𝜕𝑟
𝜇
𝜕𝑧 𝜕𝑧
S
S

(1.39)

よって式(1.25)の弱形式化が行えた．
以上より式(1.1)を弱形式化して整理すると以下のように表せる．
{𝑁}𝐴𝐴 𝜕{𝑁}
𝑟 𝜕{𝑁} 𝜕𝐴d𝑟d𝑧 + ∬d𝑟d𝑧 + ∬
d𝑟d𝑧
𝜇
𝜕𝑟
𝜇
𝜕𝑟
𝜕𝑟
S
S
S 𝜇𝑟

∬

{𝑁} 𝜕𝐴𝑟 𝜕{𝑁} 𝜕𝐴+∬d𝑟d𝑧 + ∬d𝑟d𝑧 + ∬{𝑁}𝜎(𝐴̇ + ∇𝜙)d𝑟d𝑧 = 0
𝜇
𝜕𝑟
𝜇
𝜕𝑧 𝜕𝑧
S
S
S

1.7

(1.40)

ガラーキン法ガラーキン法とは、有限要素法を行う際に微分方程式の近似解を導出するための方法である．また同様に微分方程式の近似解を導出する方法の一つである重み付き残差法を発展した方法となるため、近似解、重み付き残差法、ガラーキン法の順番で説明する．
1.7.1

近似解

近似解は解析解(厳密解)であるための条件のうちいくつかを弱くした条件を満足するものである．解析解であるための必要条件は

必要な階数だけ微分が可能である．



境界条件を満足する．



支配方程式をいたるところで満足する．

の 3 つである．
近似解の原型について未知関数を𝑢(𝑥)、力を𝑓(𝑥)とした微分方程式を
d2 𝑢(𝑥)+ 𝑓(𝑥) = 0
d𝑥 2(1.41)とし、境界条件を斉次ディレクレ境界条件
𝑢(𝑥 = 0) = 0, 𝑢(𝑥 = 𝑙) = 0とした問題の近似解を導出することで説明する．

8(1.42)

まず、未知関数𝑢(𝑥)の近似関数を有限個の異なる既知関数𝑔𝑖 (𝑥)の線形和で表す．
𝑁

𝑢̃(𝑥) = ∑ 𝑎𝑖 𝑔𝑖 (𝑥)(1.43)𝑖=1

ここで、𝑔𝑖 (𝑥)( 𝑖 = 1、2、…、N )は、解𝑢(𝑥)と同じ斉次境界条件
𝑔𝑖 (0) = 𝑔𝑖 (𝑙) = 0(1.44)

を満たす N 個の既知のなめらかな関数である．なお、ここでは𝑔𝑖 (𝑥)のことを基底関数と呼ぶ．
1.7.2

重み付き残差法

方程式の誤差を重み付きの意味で0にしようとする方法を重み付き残差法という．式(1.19)の𝑢̃(𝑥)が厳密解であると仮定すると、このとき任意の関数𝑢(𝑥)に対して、
1

∫ (
0
1

d2 𝑢̃(𝑥)+ 𝑓(𝑥)) 𝑣(𝑥)d𝑥 = 0
d𝑥 2

𝑁∫ (∑ 𝑎𝑖0𝑖=1

d2 𝑔𝑖 (𝑥)+ 𝑓(𝑥)) 𝑣(𝑥)d𝑥 = 0
d𝑥 2(1.45)(1.46)という式が成り立つ．できるだけ多くの𝑣(𝑥)について上式を満足するとより精度の良い近似解といえる．また、このような意味で関数𝑣(𝑥)のことを試験関数または重み関数と呼ぶ．
1.7.3

ガラーキン法

重み付き残差法の式(1.45)、(1.46)において𝑣𝑖 (𝑥) = 𝑔𝑖 (𝑥)とする方法をガラーキン法という．
解𝑢(𝑥)を近似するために用いた基底関数𝑔𝑖 (𝑥)そのものを試験関数とする重み付き残差法である．

1.8

グリーンの法則

局面 S の周囲を C とするとき以下の式が成り立つ
𝜕𝜙 𝜕𝜓
∬ ( − ) d𝑥d𝑦 = ∮(𝜓d𝑥 + 𝜙d𝑦)
𝜕𝑦
S 𝜕𝑥
C
1.8.1

(1.47)

ガウスの発散定理

グリーンの法則証明の際にガウスの発散定理を用いるので以下に証明する．3 次元ユーク
リッド空間の有界な領域を V、その境界面を S とする．この閉領域で定義されたベクトル関数
𝐴 (𝑥、𝑦、𝑧)、単位法線ベクトル𝒏について以下の式が成り立つ．
∭ ∇ ⋅ 𝑨d𝑉 = ∬𝑨 ⋅ 𝑛d𝑆v

S(1.48)

これは（領域全体での増減）＝（領域表面で出たり入ったりした量の差）を示している．

9

また、成分表示すると以下で表せる．
𝜕𝐴1 𝜕𝐴2 𝜕𝐴3∭(++) d𝑥d𝑦𝑑𝑧 = ∬(𝐴1 𝑐𝑜𝑠𝛼 + 𝐴2 𝑐𝑜𝑠𝛽 + 𝐴3 cosγ)d𝑆
𝜕𝑦
𝜕𝑧
v 𝜕𝑥
S

(1.49)𝑨 = (𝐴1 、𝐴2 、𝐴3 ) = (𝐴1 (𝑥、y、z)、𝐴2 (𝑥、y、z)、𝐴3 (𝑥、y、z))(1.50)𝒏 = (cos𝛼、cos𝛽、cos𝛾)(1.51)

ここで、とする．
[証明]𝑥軸に平行な直線を引き、閉曲面 S との交点を原点に近い側から𝑀1、𝑀2とする．このとき V の𝑦𝑧平面への射影をSyzとすると、体積分の𝑥軸方向成分は以下で表せる．
𝜕𝐴1
𝜕𝐴1
d𝑉 = ∬ (∫) d𝑥1 = ∬ [𝐴1 (𝑀2 ) − 𝐴1 (𝑀1)]𝑆𝑦𝑧
𝜕𝑥1
v 𝜕𝑥
S𝑦𝑧
S𝑦𝑧

∭

(1.52)

ここで、面積素d𝑆yz は、点𝑀2における曲面 S の面積素d𝑆(𝑀2)の射影として以下のように表現できる．
d𝑆yz = d𝑆(𝑀2 )cos𝛽𝑀2(1.53)

𝜕𝐴1
d𝑉 = ∬𝐴1 cos𝛽d𝑆
v 𝜕𝑥
S(1.54)

これより、∭同様に𝐴2、𝐴3軸方向成分でも成り立つので、これらの式を足すことで定理が成り立つ，[証明終わり]1.8.2

グリーンの法則の証明[証明]式(1.48)において、𝐴 = [𝜙(𝑥、𝑦)、 − 𝜓(𝑥、𝑦)、0](1.55)とし、積分範囲は底面が𝑥𝑦平面内にある領域 D(この周囲を C とする)で高さが 1 の鉛直に立つ柱とする．この柱の内部を V、表面を S=D(底面)+D’(上面)+S’(側面)とすると、式(1.48)は、
𝜕𝜙 𝜕𝜓∭ ∇ ⋅ 𝐴d𝑉 = ∭ ( − ) d𝑥d𝑦d𝑧
𝜕𝑦
v
v 𝜕𝑥
1
𝜕𝜙 𝜕𝜓
= ∬ ( − ) d𝑥d𝑦 ∫ d𝑧
𝜕𝑦
S 𝜕𝑥
0

(1.56)(1.57)

一方、式(1.48)の右辺は上面と底面に関する積分は打ち消しあう為、側面 S’についてのみ
考える．このとき、側面では、曲線 C の微分長さをd𝑠とするとd𝑆 = d𝑧d𝑠とおけるので、側面上の法線単位ベクトルを、𝒏 = (cosα、cosα、0)とすると

10(1.58)

cosαd𝑠 = dy, sinαd𝑠 = −d𝑥より∬ 𝐴 ⋅ 𝑛d𝑆 = ∬(𝜙cos𝛼 − 𝜓sin𝛼)d𝑧d𝑠
SS1= ∮ ∫ (𝜙cos𝛼 − 𝜓sin𝛼)d𝑧d𝑠
C 0= ∮(𝜙cos𝛼 − 𝜓sin𝛼)d𝑠C= ∮(𝜙d𝑦 + 𝜓d𝑥)C

(1.59)したがって
𝜕𝜙 𝜕𝜓
∬ ( − ) d𝑥d𝑦 = ∮(𝜓d𝑥 + 𝜙d𝑦)
𝜕𝑦
D 𝜕𝑥
C

(1.60)となり、D と S を改めると定理になる．
[証明終わり]1.9

磁束密度 B の計算

外部磁場がない半径 R の真っ直ぐな超伝導円柱において、長さ方向に交流電流を流す場合を考える．このとき電流により発生する自己磁場𝐻I は、交流電流の大きさを𝐼(𝑡)とすると、
𝐻I =

𝐼
2𝜋𝑅

(1.61)となる．次にパラメーター𝛼cと𝛾を用いたピン力𝐹p (𝐵̂)の関数形のモデルである Irie-Yamafujiモデル
𝐹p (𝐵̂) = 𝛼𝑐 𝐵̂𝛾(1.62)

を仮定し、方位角方向の磁束密度を𝐵、その絶対値を𝐵̂、符号因子を𝛿とすると、準静的過程における力の釣り合いの式は、−

𝐵̂ d⋅ (𝑟𝐵̂) = 𝛿𝛼𝑐 𝐵̂𝛾
𝜇0 𝑟 d𝑟(1.63)と表せる．これを解くと磁束分布は、
2−𝛾
̂𝐼 )2−𝛾 +δ(𝑟𝐵̂)= 𝛿𝑅 (𝑅𝜇0 𝐻

2−𝛾𝛼 𝜇 (𝑅 3−𝛾 − 𝑟 3−𝛾 )3−γ 𝑐 0̂𝐼 = |𝐻𝐼 |であり、𝛿R は超伝導体表面𝑟 = 𝑅におけるδの値である．
となる．ここで𝐻

11(1.64)

1.10

PHOTO-Series

PHOTO-Series は、株式会社フォトンが開発した解析ソフトウェアであり、電磁現象を利用した製品や部品などを、コンピュータ上でモデル化し有限要素法を用いて解析できる．このソフトウェアはモータや発電機などの分野や、超伝導のようなハイテク分野でよく利用されている．解析の流れとしては、まずモデルを作製し、次に物性値などを与えてから解析する．そして解析結果はコンタ図やベクトル図などの画像で表すことができる．また、デモ版（評価体験版）であれば、フリーライセンスなので誰でも気軽に扱うことができる．さらに、デモ版と製品版との違いはモデルを作製するときの要素数に制限がある程度で、解析できる内容は製品
版と何も変わりはないため、デモ版でも十分に解析を行うことができる．

1.11 本研究の目的
超伝導体について理解する時において、超伝導体内部の電磁現象の様子は対称性の高い場
合を除いてなかなか想像しづらいという問題がある．しかし、電流を流した時や磁石を近づけた時の超伝導体内部の電磁現象の様子を手計算で算出するのは複雑で非常に困難である．そのため、その問題を解決するために有限要素法を使うという方法がある．有限要素法とは解析的に解くことが難しい偏微分方程式の近似解を得るために有効な方法の一つである．したがって、
有限要素法を使って解析し超伝導体内部の電磁現象の様子を可視化することができれば、超伝導体への理解がより深まると考えられる．
本研究では手軽に解析が行える PHOTO-Series という解析ソフトウェアを用いて超伝導体を使ったモデルを作製し、超伝導体の電磁現象を有限要素法によって数値解析を行うことを目的とする．

12

第 2 章 解析方法
本研究では、超伝導体をさまざまなモデルで考えて解析を行う．ここで解析を行う理由は、酸化物超伝導体が発見され、それにはバルク超伝導体という利用方法があり、そのバルク超伝導体に磁石を近づけると磁石は浮上するが、その時の磁束密度などの様子を見たいからである．
しかし、解析モデルの全体図をそのまま解析モデルにしようとすると要素数が多くなり、また解析時間も長くかかってしまうことから、解析モデルは対称性を利用して全体図から一部切り出したものを設定した．

2.1

円筒形超伝導体に電流を軸方向に流す場合

まず 1 つ目の解析は、円筒形の超伝導体に電流を流した時の磁束密度分布を解析する．解析モデルの寸法は、円筒超伝導体の内径は 1 mm で高さは 0.5 mm に設定した．また、周りの空気の範囲は半径 3 mm、高さ 0.5 mm であり、全体図から切り出した角度は 5°で考えた．
図 2.1 の左図にモデルの全体図、右図に解析モデルを表した図をそれぞれ示す．さらに図 2.2に円筒形超伝導体の解析モデルの寸法を表した図と、図 2.3 に実際に PHOTO-Series を用いて作製したモデルを示す．次に実際に電流𝐼を流した時に磁束密度𝐵がどのように分布するかを表す図を図 2.4 に示す．また超伝導体に下から上に垂直に流れる交流電流𝐼1 = 𝐼m sin𝜔𝑡の最大電流値𝐼m を𝐼m = 3.0 × 102 A に設定した．その時のsin𝜔𝑡の波形を図 2.5 に示す．

空気
超伝導体

図 2.1: 円筒形超伝導体の全体図と解析モデル

13

3 mm

0.5 mm1mm

0.3 mm5°

図 2.2: 図 2.1 の解析モデルの寸法

図 2.3: PHOTO-Series を用いて作製した円筒超伝導体の解析モデル

14

I

B

図 2.4: 電流𝑰を流した時の磁束密度𝑩の分布の様子

図 2.5: 超伝導体に流れる交流電流の𝐬𝐢𝐧𝝎𝒕の波形

15

2.2

コイルに電流を流す場合

2 つ目の解析は、超伝導体の上にコイルがある状態をモデルとし、コイルに電流を流した
場合について行う．解析モデルの寸法は、超伝導体を半径 2 mm で高さ 1 mm、コイルは外径2 mm、内径 1 mm、高さ 1 mm であり、コイルと超伝導体の距離は 1 mm に設定した．さらに、周りの空気は半径 4 mm、高さ 5 mm で考え、また全体図から切り出した角度は 10°にした．図 2.6 の左図にモデルの全体図を、右図に解析モデルを表した図をそれぞれ示す．また図
2.7 に 解 析 モ デ ル の 寸 法 を 表 し た 図 を 示 す ． コ イ ル で 時 計 回 り に 流 れ る 直 流 電 流 𝐼2 を
𝐼2 = 3.0 × 104 A に設定し、コイルを超伝導体に近づけた時と近づけなかった時の超伝導体内の磁束密度の様子を解析した．また実際に PHOTO-Series を用いて作製した解析モデルを図2.8 に示す．さらにコイルを超伝導体に近づけた時の、コイルと超伝導体の時間毎のお互いの距離を図 2.9 に示す．このモデルは実際の浮上実験のスケールとは異なるが、本質的には問題
ない．

空気

コイル

超伝導体

全体図

解析モデル

図 2.6: 超伝導体の上にコイルがある状態のモデル

16

4mm

1mm
5mm
1mm

1mm

1mm

10°

2mm

図 2.7: 図 2.6 の解析モデルの寸法

図 2.8: PHOTO-Series を用いて作製した解析モデル17図 2.9: 時間毎のコイルと超伝導体の距離18

第 3 章 結果と考察
3.1

2.1 のモデルにおける磁束密度 B の分布

2.1 のモデルにおいて、超伝導体に交流電流を流した場合の𝑡 = 1.0 s の円筒超伝導体の磁束密度ベクトルを PHOTO-Series を用いて表した図を図 3.1 に示す．図 1 において緑色の線で囲った部分が超伝導体である．超伝導体の伝導率を計算するのに設定したパラメーターは、
𝐸0 = 1.0 × 10−4 V/m、𝐽c = 1.0 × 1010 A/m2、𝑛 = 20、𝜎max = 𝜎init × 100 1/Ωm とした．
図 3.1 を見ると、超伝導体の外側から磁束が侵入していることが分かる．また磁束密度の大きさは超伝導体の表面付近が他の領域と比べて大きくなっている．これは流した電流が超伝導体の表面に流れるからである．

図 3.1: PHOTO-Series を用いて作製したモデルで円筒超伝導体の磁束密度ベクトルの様子

19

次に超伝導体に交流電流を流した時の時間毎の超伝導体内の磁束密度を図 3.2 に示す．図3.2 を見ると時間が経つにつれて超伝導体内に侵入した磁束の密度が増加し、また磁束密度が増加する超伝導体内側からの距離の範囲が増えている．これによって時間が経つにつれてより多くの磁束が超伝導体に深く侵入していくのが確認できる．

磁束密度 B [T]3

t = 1.0 [s]t = 2.0 [s]t = 3.0 [s]

2

1

0
0

0.00010.00020.0003円筒超伝導体内側からの距離 [m]図 3.2: 円筒超伝導体内の時間毎の磁束密度20次に 2.1 のモデルにおいて、𝑡 = 3.0 s の時の𝑛値を 20、200 と変えた場合の超伝導体内の磁束密度を図 3.3 に示す．図 3.3 を見ると、𝑛値が 20 の時より 200 の時の方が磁束の超伝導体に侵入する範囲は狭いのがわかる．これは𝑛値が大きい場合、電流値が少し下がった時に発生する電圧が𝑛値の小さいときよりも急速に小さくなるため、結果的に𝑛値が 20 よりも 200 のほうが、磁束密度が減少し始める距離が長くなるからである．

n = 20n = 200

磁束密度 B [T]

2

1

0
0

0.00010.0002円筒超伝導体内側からの距離 [m]図 3.3: 𝒏 = 𝟐𝟎, 𝟐𝟎𝟎の時の磁束密度21さらに 2.1 のモデルの𝑡 = 3.0 s の場合において、σmax⁄𝜎init が 10、100、1000 の時の超伝導体内の磁束密度を図 3.4 に示す．図 3.4 を見ると、σmax⁄𝜎init が 10、100、1000 と増加するにつれて、磁束が超伝導体内部に侵入する範囲が狭くなる．これは図 1.2 よりσmax⁄𝜎init がより大きいほうが近似ラインと I-V 特性カーブのずれが小さくなり、より理論値に近づいていくからだと考えられる．

σmax / σinit = 10
σmax / σinit = 100
σmax / σinit = 1000
磁束密度B [T]

2

1

0
0

0.00010.0002円筒超伝導体内側からの距離 [m]図 3.4: 𝝈𝐦𝐚𝐱 ⁄𝝈𝐢𝐧𝐢𝐭 = 𝟏𝟎, 𝟏𝟎𝟎, 𝟏𝟎𝟎𝟎の時の磁束密度

22

3.2

2.2 のモデルにおける磁束密度 B の分布

2.2 のモデルにおいてコイルに直流電流を流し、コイルを任意の速度で超伝導体に近づけた時の𝑡 = 5.0 s の磁束密度ベクトルを PHOTO-Series を用いて表した図を図 3.5 に示す．図3.5 において赤い線で囲った部分がコイルで、茶色の線で囲った部分が超伝導体を表している．
こ の 時 設 定 し た パ ラ メー タ ー は 、 𝐸0 = 1.0 × 10−4 V/m 、 𝐽c = 1.0 × 1010 A/m2 、 𝑛 = 20 、
𝜎max = 𝜎init × 100 1/Ωm とした．
図 3.5 を見ると、コイルに電流を流したことによって発生した磁束が、超伝導体に侵入している様子が分かる．またコイルに発生する磁束は、コイルの外側よりも内側の方が大きくなっている．これはコイルが円盤状であることにより、内側の表面積より外側の表面積が大きくなることから、結果的にコイルの電流密度が、外側より内側の方が大きくなるからである．

図 3.5: PHOTO-Series を用いて作製したモデルで超伝導体にコイルを近づけた時の磁束密度ベクトルの様子

23さらにコイルを超伝導体に近づけなかった時の超伝導体内の時間毎の磁束密度を図 3.6 に、コイルを超伝導体に近づけた時の超伝導体内の時間毎の磁束密度を図 3.7 にそれぞれ示す．
図 3.6 と図 3.7 を比較すると、コイルを超伝導体に近づけた方が近づけなかった方に比べて超伝導体内の磁束密度が大きい．また𝑡 = 5.0 s に近づいていくほど磁束密度は大きくなるが、その増え方はコイルを近づけなかった方よりも、近づけた方が大きくなる．

磁束密度B [T]3t = 1.0 [s]t = 2.0 [s]t = 3.0 [s]t = 4.0 [s]t = 5.0 [s]

2

1

0
0

0.001
超伝導体中心からの距離 [m]0.002

図 3.6: コイルを超伝導体に近づけなかった時の超伝導体内の磁束密度

24磁束密度B [T]3

t = 1.0 [s]t = 2.0 [s]t = 3.0 [s]t = 4.0 [s]t = 5.0 [s]

2

1

0
0

0.001
超伝導体中心からの距離 [m]0.002

図 3.7: コイルを超伝導体に近づけた時の超伝導体内の磁束密度25第 4 章 まとめ
超伝導体内の電磁現象について有限要素法を用いた PHOTO-Series という解析ソフトウ
ェアを使用して、2 種類の超伝導体のモデルを作製した．1 つ目のモデルは円筒超伝導体のモデルで、それに交流電流を流しその時の磁束密度分布の様子を解析した．さらに、超伝導体の電気導電率を計算するプログラムのパラメーターの𝑛値と𝜎max⁄𝜎init を変更した場合の磁束密度も解析した．2 つ目のモデルはコイルと超伝導体を使用したモデルで、コイルに直流電流を流し、またそのコイルを超伝導体に近づけた場合と近づけなかった場合の、磁束密度分布の様子を解析した．
今回の解析において、まず 1 つ目の円筒超伝導体のモデルでは、超伝導体に交流電流を流して発生した磁束は、時間が経つにつれて徐々に超伝導体内に侵入していく様子が確認できた．
また同じモデルにおいて、超伝導体の電気導電率を計算するために必要なパラメーターである
𝑛値を変更した場合においては、𝑛値を大きくした方が磁束密度の発生する円筒超伝導体内側からの距離が長くなることを確認した．さらに同じ超伝導体の導電率を求めるために必要なパラメーターである𝜎max⁄𝜎init を変更した場合は、値を大きくするにつれて𝑛値の時と同じように磁束密度が発生する超伝導体内側からの距離が長くなった．これは図 1.2 より近似するカーブが超伝導体の I-V 特性カーブに近づくことから、𝜎max ⁄𝜎init を大きくしていくにつれて理論値に近づいていくと考えられる．次に 2 つ目のモデルであるコイルと超伝導体を用いたものでは、
PHOTO-Series においてコイルに直流電流を流した時に発生した磁束が、コイルを近づけた場合と近づけなかった場合に関わらず、時間が経つにつれて超伝導体内に侵入していく様子が確認できた．ここでコイルを近づけなかった場合でも時間が経つにつれて超伝導体内の磁束密度が増加した理由については、今回の研究では見つけられなかったので、このことは今後の課題になると考えられる．またコイルを超伝導体に近づけた場合と近づけなかった場合での超伝導
体内の磁束密度を比較すると、コイルを近づけた場合の方が超伝導体内の磁束密度が大きく増加していることを確認できた．このことからコイルを超伝導体に近づけた方が多くの磁束が超伝導体内に侵入していると考えられる．
以上のことから、超伝導体の電磁現象を PHOTO-Series を用いてモデルを作製し、解析を行い磁束密度分布の様子を可視化することができた．次にこの研究の今後の展望としては、今
回は PHOTO-Series という解析ソフトウェアを用いたが、有限要素法を用いて計算する他の解析ソフトでも同じ解析をして、どちらがより優れているかを検証することや、今回は割合単純なモデルで解析を行ったが、次はより複雑なモデルで行い解析を現実的な状況で考えることなどが挙げられる．

26

謝辞
本研究を行うにあたり、小田部荘司教授に多大なるご指導、助言を頂き、深く感謝いたします．また、様々な助言やご指導、ご協力をして頂いた松下照男教授、木内勝准教授に深く感謝いたします．そして、研究がうまく進まなかった際に、研究の解析を手伝っていただくなど、様々なご指導を頂いた小松伸二郎さん、和田純さんをはじめ研究室内でお世話になりました小
田部研究室、木内研究室所属の皆様にも深く感謝いたします．

27

参考文献
[1] 本庄卓一、「超電導ケーブルに生じる交流損失の数値解析」、JMAG Users Conference 2001
[2] 筒井喜平、「有限要素法による磁場解析」、http://computation.cside.com/mag-index.html[3] 株式会社フォトンホームページ、http://www.photon-cae.co.jp/product/series/index.html[4]「有限要素法特論」、http://www.sml.k.u-tokyo.ac.jp/members/nabe/lecture2003/01.pdf[5] 松下照男、「—新電磁気学-電気・磁気学の新しい体系の確立—」、コロナ社、1994 年

28