Pythonでお手軽に3Dトポロジー最適化してみよう

トポロジー最適化試してみたいよね。でも・・・

「トポロジー最適化、興味はあるけど、なんだか気軽に試せない…」
そんなふうに感じたことはありませんか？私自身、最初に取り組もうとしたとき、いくつかの「見えない壁」にぶつかりました。

たとえば：

準備が大変そう：有限要素法（FEM）や構造解析の知識が必要とされるイメージが強く、手を出しづらい。

ライブラリやツールがバラバラ：解析と最適化を別々に組み合わせないといけないケースが多く、環境構築で疲れてしまう。

そもそもどうやって形を変えるの？ という根本的な疑問に対する情報が少ない。

このように、ハードルが高い印象があるのは事実です。

でも意外と簡単な一面も

ところが、一歩踏み込んでみると「意外とシンプルじゃないか？」と感じる瞬間があります。
なぜなら、**「構造解析さえ解ければ、あとは最適化するだけ」**だからです。

構造解析の解法（FEMなど）さえ手元にあれば、

・各要素の剛性をどう変えるか？
・形状をどう更新するか？
・制約条件（体積制限など）をどう満たすか？

といった問題は、ルールに沿って繰り返し計算していくだけの仕組みに落とし込むことができます。
近年では Python ライブラリ（例：scikit-fem, PETSc, DL4TO, pytopo3d など）も充実しており、オープンソースで十分に楽しめる状況と思います。

試してみよう

とはいえ、ゼロから自分で最適化アルゴリズムを実装するのは、やはり骨が折れます。
そこで今回は、Python製のトポロジー最適化ライブラリ scikit-topt を使って、
できるだけ簡単に「最適な構造」を求めてみようと思います。

このライブラリは、最低限の記述でトポロジー最適化を回すことができ、
中身もシンプルなので、自分で改造したい人」や「理論を学びたい人」にもおすすめです。
体感できれば、きっと理解も早いはず。

なお、計算の理屈については多くの人の記事が見つかると思うので、大部分は省略します。
https://qiita.com/fujitagodai4/items/7cad31cc488bbb51f895

ただし一点、複数荷重条件のトポロジー最適化について書いておきます。

複数荷重条件のトポロジー最適化とは？

トポロジー最適化では、一般的に「荷重をかけたときに変形しにくい（≒剛性が高い）」構造を見つけることを目的とします。
これは、コンプライアンス（変形エネルギー）を最小化する問題として定式化されます。

単一の荷重条件では、以下のような最適化問題になります：

\begin{aligned}
\min_{\rho} \quad & C = \mathbf{f}^\top \mathbf{u} = \sum_{e} \rho_e^p \mathbf{u}_e^\top \mathbf{K}_e \mathbf{u}_e \\
\text{subject to} \quad & \mathbf{K}(\rho) \mathbf{u} = \mathbf{f} \\
& \sum_e \rho_e v_e \leq V_{\max}, \quad 0 \leq \rho_e \leq 1
\end{aligned}

ここで：

• $\rho_e$：各要素の密度（設計変数）
• $\mathbf{K}_e$：要素剛性マトリクス
• $\mathbf{u}_e$：要素変位ベクトル
• $p$：SIMP法のペナルティ係数
• $V_{\max}$：許容体積

複数荷重条件がある場合は？

現実の構造物は、複数の異なる荷重ケースにさらされることが多いです。たとえば：

• 自動車部品なら、前後・上下の揺れに対応する必要がある
• 橋なら、車両の走行する場所が時間変化する

こうした場合、各荷重条件ごとにFEMを解き、それぞれの変形エネルギー（コンプライアンス）を合計したものを最小化します：

\begin{aligned}
\min_{\rho} \quad & \sum_{i=1}^{N} w_i \, \mathbf{f}_i^\top \mathbf{u}_i = \sum_{i=1}^{N} w_i \sum_e \rho_e^p \mathbf{u}_{i,e}^\top \mathbf{K}_e \mathbf{u}_{i,e} \\
\text{subject to} \quad & \mathbf{K}(\rho) \mathbf{u}_i = \mathbf{f}_i \quad \text{for } i = 1, \dots, N \\
& \sum_e \rho_e v_e \leq V_{\max}, \quad 0 \leq \rho_e \leq 1
\end{aligned}

ここで：

• $N$：荷重ケースの数
• $\mathbf{f}_i$, $\mathbf{u}_i$：各荷重ケースの外力・変位ベクトル
• $w_i$：各荷重ケースの重み（通常はすべて1でもよい）

トポロジー最適化のための環境構築

本記事では、scikit-toptというライブラリを用います。
https://github.com/kevin-tofu/scikit-topt

インストールには、Python環境を用意してpipするだけです。

pip install scikit-topt

実験、体感してみよう

問題を設定する

まず、扱う領域をメッシュとして作成します。
ここでは、x, y, zの長さがそれぞれ8, 8, 1の直方体を定義し、Hexahedral1次要素でメッシュを切ります。

import skfem
import sktopt

x_len = 8.0
y_len = 8.0
z_len = 1.0
mesh_size = 0.2
mesh = sktopt.mesh.toy_problem.create_box_hex(
    x_len, y_len, z_len, mesh_size
)
e = skfem.ElementVector(skfem.ElementHex1())
basis = skfem.Basis(mesh, e, intorder=2)

次に、ディリクレ境界条件を課す領域を設定します。

dirichlet_nodes = sktopt.mesh.utils.get_nodes_indices_in_range(
    basis.mesh, (0.0, 0.05), (0.0, y_len), (0.0, z_len)
)
dirichlet_dir = "all"

そして、力を課すノードを設定します。
ここでは、二つの領域にそれぞれY軸方向に100の力を加えています。
F_nodes_0にはマイナス方向に、F_nodes_1にはプラス方向に力を加えています。
このとき、物体にすべての力を一度に加えるのではなく、それぞれの荷重条件ごとにFEM（有限要素法）を解くことになります。
つまり、先に述べたように、異なる方向から個別に力がかかる場合でも、それぞれに対して十分に強い構造を作ろうと、最適化アルゴリズムが働きます。

F_nodes_0 = sktopt.mesh.utils.get_nodes_indices_in_range(
    basis.mesh, (x_len, x_len), (y_len, y_len), (0, z_len)
)
F_nodes_1 = sktopt.mesh.utils.get_nodes_indices_in_range(
    basis.mesh, (x_len, x_len), (0, 0), (0, z_len)
)
F_nodes = [F_nodes_0, F_nodes_1]
F_dir = ["u^2", "u^2"]
F = [-100, 100]

設計領域を設定します。
単純に全領域の密度を最適化対象にしています。

design_elements = sktopt.mesh.utils.get_elements_in_box(
    mesh,
    (0.0, x_len), (0.0, y_len), (0.0, z_len)
)

タスクをオブジェクトとして生成します。

E0 = 210e9
nu = 0.3
mytask = sktopt.mesh.task.TaskConfig.from_nodes(
    E0,
    nu,
    basis,
    dirichlet_nodes,
    dirichlet_dir,
    F_nodes,
    F_dir,
    F,
    design_elements
)

作成した問題を図示すると以下のようになります。

実験してみよう

さて、実験してみましょう。

cfg = sktopt.core.optimizers.LogMOC_Config(
    vol_frac=0.6,
    max_iters=40,
    record_times=40,
    export_img=True
)
optimizer = sktopt.core.LogMOC_Optimizer(cfg, tsk)
optimizer.parameterize()
optimizer.optimize()

さて、結果ですが以下のようになりました。
そこそこ、いい感じですね！

最適化に利用する様々なパラメータがあるのですが、説明は省略します。
また機会あれば・・・