1 はじめに
想定読者
- 電磁界解析ソフトを利用しているが，その中身を知らない人
- 電磁気学をある程度理解しており，シミュレーションに興味がある人
- 私の友人O君
方針
- 最も基本的な2次元静磁場の電磁界解析の実装を理解できるようになる．特に今回は3角形メッシュを用いる．
- 環境構築が簡単なPythonで実装する．Pythonは機械学習をはじめとする科学分野で広く使用される．
- マスクウェル方程式から実装方法を説明する．ただし，詳細な式変換や物理的な意味の解説は別の資料に譲る．あくまでも実装に必要な部分の詳細に限る．
- この後の発展内容についての参考資料を示す．
背景
電磁界解析は現在の生産技術を支える礎となっています．モーターやインダクタ等の電気機器を解析することによりトルクや発熱、振動がわかります．電気機器の設計を行っている方であればJMAG等の商用ソフトを使用してそれらの特性を解析していることでしょう．
さて、それらのソフトで実行されている電磁界解析について設計者は理解できているでしょうか．
私の経験から言わせていただきますと、その仕組みを知る人は少数です．
電磁界解析のとっかかりは難しく，理解をしようとしてはみても挫折した人は多いのではないでしょうか．ググってもヨクワカラナイ式の海に飲まれます．さらに悪いことに、簡単なソースコードさえ見つけられない．
今回の記事は基本である有限要素法の支配方程式を導出し実装します．マクスウェル方程式から始まり、磁場を可視化できるまでを案内します．
なお、詳細な定式化や物理的な意味は既存の資料に譲り、ここでは読者の手元で電磁界解析を動かせるようになることに注力します．
ぜひ，手元にソースコードをダウンロードして，その仕組みを追ってみてください．
コードの取得
ファイルを取得したいフォルダで下記コマンドを実行してください．
git clone https://github.com/kitsune8848/001_FEM.git
ファイル取得後は001_FEMフォルダに入っているmain.pyを実行することで電磁界解析が実行されます．
2 解析対象
解析対象のモデル
今回はとても簡単なモデルを用意しました．
鉄の周りにコイルがあるシンプルなモデルです．コイルに電流が流れることで磁束が発生します．
Z方向は同一の形状であり影響が無視できるくらい長いと仮定すると2次元解析で扱えます．
発生する時速の様子をなんとなく想像できますでしょうか．磁束は透磁率の高い鉄を集中して通るはずですよね？それをシミュレーションにより確認します．
解析するためのメッシュモデル
PCで対象を解析するためには有限要素法と呼ばれる解析対象を微表な要素に分割する方法を用います．この方法は電磁気に関わらず橋の設計や熱のシミュレーションで広く用いられています．各要素もしくは節点で特性を求めます．
解析モデルはフォルダ内のinput_model.vtkです．可視化には無料ソフトのparaviewを使用します．なお，解析モデルの対称性より4分の1モデルで十分です.
なお，このモデルは要素数2732, 節点数1442です．
解析結果
この後解説する電磁界解析の結果を先に紹介します．Pythonを実行しますとフォルダ内にinput_model.vtkが作成されますので，paraviewで可視化してみてください．次の図のように磁束密度が得られています．以下のような特徴がちゃんと得られていますね．
- 透磁率の高い鉄に磁束が集中している
- 電流が流れるコイルの周り磁束ができている
3 定式化
定式化のゴール
まず，今回の定式化のゴールを示しておきましょう．各節点上のベクトルポンテシャル$A_i$のベクトル$\mathbf{a}$を求めることです．今回の場合だと節点が1442あるので1442次元の行列を扱うことになります．
\mathbf{a} =
\begin{bmatrix}
A_1 \\
A_2 \\
\dots \\
A_n
\end{bmatrix} \tag{1}
ベクトル$\mathbf{a}$は剛性行列と呼ばれる$\mathbf{M}$と右辺のベクトル$\mathbf{b}$から次のような式が成り立ちます．
$$
\mathbf{M} \mathbf{a} = \mathbf{b} \tag{2}
$$
マクスウェル方程式を変形することにより，合成行列と右辺ベクトルを導出します．
これらが，わかれば$\mathbf{M}$の逆行列を求めることにより，$\mathbf{a}$は分かりますね．(厳密には有限要素法では逆行列を直接求めることはせず，cg法を用います．)
ベクトルポテンシャルがわかると磁束密度を求められます．磁束密度がわかればマスクウェル応力から電磁力がわかります．さらに，振動や熱解析も可能になります．
2次元静磁場における支配方程式の導出
みなさまご存知マスクウェル方程式から初めて，Pythonで実装できるレベルまで式を整理していきます．
電磁気は4つのマクスウェル方程式に支配されます．
$$
\nabla \cdot \mathbf{D} = \rho \tag{3.1}
$$
$$
\nabla \cdot \mathbf{B} = 0 \tag{3.2}
$$
$$
\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}\tag{3.3}
$$
$$
\nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}\tag{3.4}
$$
今回は簡単のため静磁場を仮定します．周波数が低い場合や誘導電場を無視できる場合は静磁場と仮定できます．トランスやモーターの解析では静磁場としてよくモデル化されます．(3.3), (3.4)は次に変形できます．
$$
\nabla \times \mathbf{E} = 0 \tag{4}
$$
$$
\nabla \times \mathbf{H} = \mathbf{J_0} \tag{5}
$$
ここで$\mathbf{J_0}$は強制電流密度ベクトルです．
以上のマスクウェル方程式をゴチャゴチャすると，2次元静磁場問題におけるベクトルポテンシャルの場の支配方程式(6)が得られます．静磁場は(6)式に従って分布が決定されます．と言ってもどんな分布になるか，わかる人はまずいないです．
$$
\nabla \cdot \left( \nu \nabla A_z \right) = -J_0\tag{6}
$$
3.2 支配方程式の離散化
静磁場では任意の座標で(6)が成り立ちます．しかし有限のメモリのPCで無限次元を扱うことは不可能ですので離散化する必要があります．そのためのメッシュです．
有限要素法では，支配方程式を離散化する際に誤差が生ます.つまり，有限の座標しか扱えないので，任意の座標において(6)を完全に満たした分布をPCで得るのは難しいということです．それではどうするか．この誤差を重み付き平均和において 0 とすることで解決できます．方法を重み付き残差法と呼ぶ.ここで重み関数 𝒘 を導入し，重み付き残差法を適用することで，次式が得られます．
$$
\int_{\Omega} w \nabla \cdot (\nu \nabla A_z) dS = - \int_{\Omega} w J_0 \ dS\tag{7}
$$
この式は２段階微分になっており，解く上で非常に扱いにくいです．微分の回数を下げ，解きやすい形にします．変形過程は省略しますが弱形式と呼ばれる式が次になります．
$$
\int_{\Omega} \nabla w \cdot (\nu \nabla A_z) dS = - \int_{\Omega} w J_0 \ dS\tag{8}
$$
ここで，今回は3角形の要素を用います．3角形の頂点3点のベクトルポテンシャルの値から３角形内の任意の座標でのベクトルポテンシャルを補完関数$N_i$を用いて求めます．ぱっと見難しそうですが，ようは3点のベクトルポテンシャルの値を通るような平面で補完しています．
$$
A_z(x, y) = \sum_{i=1}^{3} N_i(x, y) A_i\tag{9}
$$
$$
N_i(x, y) = \frac{b_i + c_i x + d_i y}{2S_a}\tag{10}
$$
\begin{aligned}
b_1 &= x_2 y_3 - x_3 y_2, &\quad c_1 &= y_2 - y_3, &\quad d_1 &= x_3 - x_2 \\
b_2 &= x_3 y_1 - x_1 y_3, &\quad c_2 &= y_3 - y_1, &\quad d_2 &= x_1 - x_3 \\
b_3 &= x_1 y_2 - x_2 y_1, &\quad c_3 &= y_1 - y_2, &\quad d_3 &= x_2 - x_1
\end{aligned}\tag{11}
ここで$S_a$は三角形要素の面積，$x_i,y_i,z_i$はその三角形要素の各頂点（節点）の座標です
さて，まだ重み$w$をどのように決めるかという問題が残っていますね．これは初学者には，もう受け入れろという他ないのですが，$w$には補完関数$N_i$を使用します．これはガラーキン法と呼ばれる有限要素解析の中では非常に重要な方法です．
離散化とガラーキン法を(8)に適用するとある要素に対して次の式が得られます．(14)を全ての要素に対して計算し全体の剛性行列$\mathbf{M}$に足し込むことで$\mathbf{M}$が得られます．
$$
\int_{\Omega} \nabla N_j \cdot (\nu \nabla \sum_{i=1}^3 N_i A_i) , dS = - \int_{\Omega} N_j J_0 \ dS \quad \text{where } j = 1, 2, 3\tag{12}
$$
(12)を行列形式にすると，
\int_\Omega \nu
\begin{bmatrix}
\nabla N_1 \cdot \nabla N_1 & \nabla N_1 \cdot \nabla N_2 & \nabla N_1 \cdot \nabla N_3 \\
& \nabla N_2 \cdot \nabla N_2 & \nabla N_2 \cdot \nabla N_3 \\
\text{sym.} & & \nabla N_3 \cdot \nabla N_3
\end{bmatrix}
\begin{bmatrix}
A_1 \\
A_2 \\
A_3
\end{bmatrix}
d S=
-J_0\int_\Omega
\begin{bmatrix}
N_1 \\
N_2 \\
N_3 \\
\end{bmatrix}
d S\tag{13}
(13)の各項を計算すると
\frac{1}{4S_a \mu_0\mu_r}
\begin{bmatrix}
c_1 c_1 + d_1 d_1 & c_1 c_2 + d_1 d_2 & c_1 c_2 + d_1 d_2\\
& c_2 c_2 + d_2 d_2 & c_2 c_3 + d_2 d_3 \\
\text{sym.} & & c_3 c_3 + d_3 d_3
\end{bmatrix}
\begin{bmatrix}
A_1 \\
A_2 \\
A_3
\end{bmatrix}
=
-
\frac{J_0 S_a}{3}
\begin{bmatrix}
1 \\
1 \\
1 \\
\end{bmatrix}\tag{14}
4 Pythonで実装
剛性行列
詳細なコードはgitにあげていますので，そちらを参照し，実行してください．行列を作成する重要な部分だけこちらに載せます．
まずは剛性行列$\mathbf{M}$を求めるコードです．剛性行列$\mathbf{M}$は次のローカルな剛性行列$\mathbf{M_L}$を要素の数だけ足し合わせて作成します．
\mathbf{M_L} =
\frac{1}{4S_a\mu_0 \mu_r}
\begin{bmatrix}
c_1 c_1 + d_1 d_1 & c_1 c_2 + d_1 d_2 & c_1 c_2 + d_1 d_2\\
& c_2 c_2 + d_2 d_2 & c_2 c_3 + d_2 d_3 \\
\text{sym.} & & c_3 c_3 + d_3 d_3
\end{bmatrix}
class Fem:
def __init__(self, mesh: Mesh):
self.mesh = mesh # メッシュ情報（節点・要素を保持）
self.current_density = 3e6 # 電流密度 [A/m^2]
self.mu = 4.0 * np.pi * 1.0e-7 # 真空の透磁率 μ₀ [H/m]
def _make_left_hand_side_matrix(self):
"""
FEMの剛性行列を生成
"""
num_nodes = len(self.mesh.nodes)
M = np.zeros((num_nodes, num_nodes))
index = [[1, 2], [2, 0], [0, 1]] # 順番指定
for element in self.mesh.elements.values():
x = [self.mesh.nodes[nid].x for nid in element.node_ids]
y = [self.mesh.nodes[nid].y for nid in element.node_ids]
ce = np.zeros(3)
de = np.zeros(3)
for i in range(3):
ce[i] = y[index[i][0]] - y[index[i][1]]
de[i] = x[index[i][1]] - x[index[i][0]]
# 要素ごとの局所剛性行列を計算
M_local = np.zeros((3, 3))
coef = 1.0 / (self.mu * element.mu_r * 4.0 * element.area)
for i in range(3):
for j in range(3):
M_local[i][j] = coef * (ce[i] * ce[j] + de[i] * de[j])
# 全体剛性行列へアセンブリ
for i in range(3):
ii = element.node_ids[i]
for j in range(3):
jj = element.node_ids[j]
M[ii, jj] += M_local[i][j]
return M
右辺ベクトル
次に右辺ベクトル$\mathbf{b}$です.こちらも同様にローカルな右辺ベクトル$\mathbf{b_L}$をy要素の数だけ足し込んで生成します．
\mathbf{b_L}
=
-
\frac{J_0 S_a}{3}
\begin{bmatrix}
1 \\
1 \\
1 \\
\end{bmatrix}
class Fem:
def __init__(self, mesh: Mesh):
self.mesh = mesh # メッシュ情報（節点・要素を保持）
self.current_density = 3e6 # 電流密度 [A/m^2]
self.mu = 4.0 * np.pi * 1.0e-7 # 真空の透磁率 μ₀ [H/m]
def _make_right_hand_side_vector(self):
"""
電流密度による右辺ベクトル b を作成
"""
b = np.zeros(len(self.mesh.nodes))
for element in self.mesh.elements.values():
if element.material == Material.COIL:
if element.cell_type == CellType.TRIANGLE:
for node_id in element.node_ids:
b[node_id] += self.current_density * element.area / 3.0
else:
raise Exception("三角形以外の要素には未対応です")
return b
境界条件
次に境界条件を適用します．詳細は省きますが，モデルの外側の境界に位置する節点のうち，磁束が外に出て行かない場合は一定のベクトルポンテシャルを与えます．よって境界条件を適用する節点のベクトルポテンシャルは求める必要がないので$\mathbf{M}, \mathbf{b}$を整理します．
class Fem:
def __init__(self, mesh: Mesh):
self.mesh = mesh # メッシュ情報（節点・要素を保持）
self.current_density = 3e6 # 電流密度 [A/m^2]
self.mu = 4.0 * np.pi * 1.0e-7 # 真空の透磁率 μ₀ [H/m]
def _apply_boundary_condition(self, A: np.ndarray, b: np.ndarray, val: float):
"""
Dirichlet境界条件を適用（ポテンシャルを固定）
"""
essential_node_ids = [
node.id for node in self.mesh.nodes.values() if node.is_essential_boundary_node
]
n = A.shape[0]
# 他の節点の式から影響を削除
for j in range(n):
for i in essential_node_ids:
b[j] -= A[j, i] * val
A[j, i] = 0.0
# 対象節点に境界値を設定（単位行列）
for i in essential_node_ids:
A[i, :] = 0.0
A[i, i] = 1.0
b[i] = val
ベクトルポンテシャルの計算
以上により$\mathbf{M}, \mathbf{b}$が分かったので，後は$\mathbf{M}$の逆行列を求めるだけです．と言っても，実はこれ難しい問題です．今回は節点約1400ですから行列の次元は1400×1400でこの逆行列を求めるのは結構大変です．さらに複雑な機器のモデルを扱うようになると数十万次元のデータを扱うことはよくあり，そうなるとマシンパワーで解決できるものではなくなります．そこで使用されるのがiccg法です．これは剛性行列が疎行列かつ正定であることを利用して効率的に解を得る方法です．今回はPythonのライブラリに頼って実装しています．
class Fem:
def __init__(self, mesh: Mesh):
self.mesh = mesh # メッシュ情報（節点・要素を保持）
self.current_density = 3e6 # 電流密度 [A/m^2]
self.mu = 4.0 * np.pi * 1.0e-7 # 真空の透磁率 μ₀ [H/m]
def analyse(self):
"""
FEM解析のメイン関数
"""
A = self._make_left_hand_side_matrix() # 剛性行列（左辺）
b = self._make_right_hand_side_vector() # 右辺ベクトル（電流源）
self._apply_boundary_condition(A, b, 0) # 境界条件を適用（Dirichlet）
x, info = MatrixSolver.iccg(A, b) # ICCG法による連立方程式の解法
さて，これで終わりと言いたいところですが，まだ磁束密度を求めていません．
磁束密度を求める
以上で得られたベクトルポンテシャルから磁束密度を求めます．
さて，磁束密度の導出は簡単なので読者に委ねます．
実装は次のようになります．
class Fem:
def __init__(self, mesh: Mesh):
self.mesh = mesh # メッシュ情報（節点・要素を保持）
self.current_density = 3e6 # 電流密度 [A/m^2]
self.mu = 4.0 * np.pi * 1.0e-7 # 真空の透磁率 μ₀ [H/m]
def _compute_flux_density(self):
"""
各要素に対して磁束密度ベクトル Bx, By を計算
"""
index = [[1, 2], [2, 0], [0, 1]]
for element in self.mesh.elements.values():
x = [self.mesh.nodes[nid].x for nid in element.node_ids]
y = [self.mesh.nodes[nid].y for nid in element.node_ids]
ce = np.zeros(3)
de = np.zeros(3)
for i in range(3):
ce[i] = y[index[i][0]] - y[index[i][1]]
de[i] = x[index[i][1]] - x[index[i][0]]
Bx = 0.0
By = 0.0
for i in range(3):
node_id = element.node_ids[i]
A_i = self.mesh.nodes[node_id].vector_potential
Bx += -A_i * de[i] / (2.0 * element.area)
By += A_i * ce[i] / (2.0 * element.area)
element.set_flux_density([Bx, By])
5 さいごに
今回2次元静磁場という最も基本的なFEMを実装しました．なんとなくわかるようなわからないような結果かもしれません．有限要素法を数時間で理解するのは至難の業です．この記事を入門として，以下の資料をあたっていただければと思います．要望があればより詳細な記事を出すかもしれません．
Q.より詳細に導出や物理的な意味を知りたい
Q. 非線形な特性を持つ磁性体を使いたい
Q. 磁石の影響を加味した解析を行いたい
A. 「電気工学の有限要素法」がおすすめです．私は初学時にこれを熟読しました．
Q. 自分でメッシュモデルを作成して解析したい
A. gmeshでググってください．
Q. 3次元の電磁界解析を行いたい
A. 「三次元有限要素法: 磁界解析技術の基礎」をあたってください．過去の私は難しくて断念しました．
Q. 電磁界解析を用いて電気機器の最適化を行いたい
A. 「電磁界解析による最適設計：トポロジー最適化の基礎から機械学習まで」おすすめです．