🌐 JP | 🇬🇧 EN | Last sync: 2025-11-16

第5章:Python実践:3Dプリンティングシミュレーション

STL解析・スライシング・熱伝導・パラメータ最適化・機械学習をコードで動かす

📚 3Dプリンティング入門シリーズ ⏱️ 読了時間: 45-50分 🎓 難易度: 中級〜上級 💻 コード例: 5 📝 演習問題: 5

学習目標

この章を完了すると、以下を説明・実践できるようになります:

基本理解(Level 1)

実践スキル(Level 2)

応用力(Level 3)

💡 この章の位置づけ

第1章から第4章では、積層造形(Additive Manufacturing, AM)の原理・材料・プロセスを体系的に学びました。本章はシリーズの集大成として、それらの知識を実際に動くPythonコードへ落とし込みます。STLの読み込みからスライシング、熱シミュレーション、パラメータ最適化、機械学習まで、5つの実践課題を順に手を動かして確認します。ここで扱うコードはすべて python3 で実行し、掲載する出力は実際の実行結果です。外部ライブラリは NumPy と scikit-learn のみに絞り、STL処理は専用ライブラリに頼らず自前で実装することで、内部の仕組みまで見通せるようにしています。

5.1 STLファイルの読み込みと解析

5.1.1 三角形メッシュとしての物体表現

STL(STereoLithography)ファイルは、物体の表面を三角形メッシュ(Triangle Mesh)の集合として表現します。各三角形は3つの頂点座標と1つの外向き法線ベクトルを持ちます。第1章で構造を学んだこの形式を、ここでは自前のコードで解析します。numpy-stlのような専用ライブラリを使わず、NumPyだけで三角形配列を組み立てることで、体積や表面積が「どの式から」得られるのかを明確にします。

解析対象として、一辺20 mmの立方体を12枚の三角形(各面2枚)として生成します。立方体は体積・表面積の理論値(それぞれ 20³ = 8000 mm³、6×20² = 2400 mm²)が分かっているため、実装の検算に最適です。

5.1.2 体積・表面積・水密性の計算

三角形メッシュから幾何量を求める式は次の3つが基本です。

コード例1: NumPyによるSTLメッシュ解析

import numpy as np

def make_cube_mesh(size=20.0):
    """一辺 size の立方体を12枚の三角形メッシュ(法線外向き)として生成"""
    s = size
    v = np.array([
        [0, 0, 0], [s, 0, 0], [s, s, 0], [0, s, 0],   # 底面 z=0
        [0, 0, s], [s, 0, s], [s, s, s], [0, s, s],   # 上面 z=s
    ], dtype=float)
    # 各面を反時計回り(外から見て CCW)で定義
    faces = [
        (0, 3, 2), (0, 2, 1),   # 底面 (-z)
        (4, 5, 6), (4, 6, 7),   # 上面 (+z)
        (0, 1, 5), (0, 5, 4),   # 前面 (-y)
        (2, 3, 7), (2, 7, 6),   # 背面 (+y)
        (1, 2, 6), (1, 6, 5),   # 右面 (+x)
        (0, 4, 7), (0, 7, 3),   # 左面 (-x)
    ]
    tris = np.array([[v[a], v[b], v[c]] for a, b, c in faces])
    return tris

def triangle_normals(tris):
    """各三角形の単位法線ベクトル (v2-v1)x(v3-v1)"""
    e1 = tris[:, 1] - tris[:, 0]
    e2 = tris[:, 2] - tris[:, 0]
    n = np.cross(e1, e2)
    lengths = np.linalg.norm(n, axis=1, keepdims=True)
    return n / lengths

def surface_area(tris):
    """三角形面積の総和 = 表面積"""
    e1 = tris[:, 1] - tris[:, 0]
    e2 = tris[:, 2] - tris[:, 0]
    return 0.5 * np.linalg.norm(np.cross(e1, e2), axis=1).sum()

def mesh_volume(tris):
    """符号付き四面体法による体積 (原点と各三角形が作る四面体の総和)"""
    v1, v2, v3 = tris[:, 0], tris[:, 1], tris[:, 2]
    signed = np.einsum('ij,ij->i', v1, np.cross(v2, v3)) / 6.0
    return abs(signed.sum())

def is_watertight(tris, tol=6):
    """各無向エッジがちょうど2枚の三角形に共有されるかを検査"""
    verts = np.round(tris.reshape(-1, 3), tol)
    uniq, inv = np.unique(verts, axis=0, return_inverse=True)
    idx = inv.reshape(-1, 3)
    edge_count = {}
    for a, b, c in idx:
        for u, w in [(a, b), (b, c), (c, a)]:
            key = (min(u, w), max(u, w))
            edge_count[key] = edge_count.get(key, 0) + 1
    counts = np.array(list(edge_count.values()))
    watertight = bool(np.all(counts == 2))
    return watertight, len(uniq), len(edge_count), counts

mesh = make_cube_mesh(20.0)
normals = triangle_normals(mesh)
area = surface_area(mesh)
vol = mesh_volume(mesh)
wt, n_vertices, n_edges, counts = is_watertight(mesh)
bbox_min = mesh.reshape(-1, 3).min(axis=0)
bbox_max = mesh.reshape(-1, 3).max(axis=0)

print(f"三角形数        : {len(mesh)}")
print(f"ユニーク頂点数  : {n_vertices}")
print(f"ユニークエッジ数: {n_edges}")
print(f"表面積          : {area:.2f} mm^2  (理論値 6*20^2 = {6*20**2})")
print(f"体積            : {vol:.2f} mm^3  (理論値 20^3 = {20**3})")
print(f"バウンディングボックス: min={bbox_min}, max={bbox_max}")
print(f"水密(watertight): {wt}  (全エッジ共有数の最小={counts.min()}, 最大={counts.max()})")
print(f"先頭3枚の単位法線:")
for i in range(3):
    print(f"  面{i}: [{normals[i,0]:+.2f} {normals[i,1]:+.2f} {normals[i,2]:+.2f}]")

実行結果:

三角形数        : 12
ユニーク頂点数  : 8
ユニークエッジ数: 18
表面積          : 2400.00 mm^2  (理論値 6*20^2 = 2400)
体積            : 8000.00 mm^3  (理論値 20^3 = 8000)
バウンディングボックス: min=[0. 0. 0.], max=[20. 20. 20.]
水密(watertight): True  (全エッジ共有数の最小=2, 最大=2)
先頭3枚の単位法線:
  面0: [+0.00 +0.00 -1.00]
  面1: [+0.00 +0.00 -1.00]
  面2: [+0.00 +0.00 +1.00]

体積8000 mm³・表面積2400 mm²が理論値と完全に一致し、実装が正しいことが確認できます。ユニークエッジ数が18・全エッジの共有数がすべて2であることは、このメッシュが水密(閉じた表面)であることを意味します。頂点8・エッジ18・面12はオイラーの多面体公式 V−E+F = 8−18+12 = 2 も満たしています。法線がすべて外向き(底面が −z、上面が +z)に揃っている点も、造形時の内外判定が正しく行える条件です。

⚠️ 水密性チェックの前提

ここでのエッジ共有チェックは、頂点座標を小数6桁で丸めて同一視しています。実際のSTLファイルでは、CADからの書き出し時に同じ位置の頂点がわずかに異なる浮動小数点値を持つ「頂点の重複」が頻繁に起こります。丸めの桁数を誤ると、本来つながっているエッジが別物と判定され、健全なメッシュを「穴あき」と誤検出します。実務では trimesh などのライブラリが許容誤差付きで頂点を統合してから判定します。ここでは仕組みを理解するために最小限の実装にとどめています。

5.2 スライシングアルゴリズムの実装

5.2.1 平面と三角形の交差

スライシング(Slicing)とは、3Dモデルを一定の高さ(レイヤー高さ)ごとに水平面で切断し、各層の輪郭を抽出する処理です(第1章参照)。その核心は「平面 z = z₀ と各三角形の交線分を求める」という幾何演算にあります。三角形の3辺のうち、平面を挟んで反対側に端点を持つ辺だけが平面と交わります。交点は端点間の線形補間で求められ、1つの三角形からはちょうど2つの交点、すなわち1本の線分が得られます。

交点 = p₀ + t·(p₁ − p₀), t = (z − z_p₀) / (z_p₁ − z_p₀)

各層でこれらの線分を集めると、その層の断面輪郭(contour)になります。ここでは断面が高さとともに変化する例として、底辺20 mm・高さ30 mmの四角錐をスライスします。四角錐の断面は正方形で、高さが上がるほど小さくなるため、輪郭抽出が正しく機能しているかを周長の減少から確認できます。

コード例2: 平面交差によるスライシング

import numpy as np

def make_pyramid_mesh(base=20.0, height=30.0):
    """底辺 base、高さ height の四角錐メッシュ(側面4枚+底面2枚)"""
    b, h = base, height
    apex = [b/2, b/2, h]
    v = np.array([
        [0, 0, 0], [b, 0, 0], [b, b, 0], [0, b, 0], apex,
    ], dtype=float)
    faces = [
        (0, 2, 1), (0, 3, 2),        # 底面
        (0, 1, 4), (1, 2, 4),        # 側面
        (2, 3, 4), (3, 0, 4),
    ]
    return np.array([[v[a], v[b], v[c]] for a, b, c in faces])

def slice_at_z(tris, z):
    """平面 z=z0 と各三角形の交線分を求め、輪郭周長を返す"""
    segments = []
    for tri in tris:
        pts = []
        for i in range(3):
            p0, p1 = tri[i], tri[(i + 1) % 3]
            z0, z1 = p0[2], p1[2]
            if (z0 - z) * (z1 - z) < 0:          # エッジが平面をまたぐ
                t = (z - z0) / (z1 - z0)
                pts.append(p0 + t * (p1 - p0))
        if len(pts) == 2:
            segments.append((pts[0], pts[1]))
    perimeter = sum(np.linalg.norm(a - b) for a, b in segments)
    return segments, perimeter

pyr = make_pyramid_mesh(base=20.0, height=30.0)
layer_height = 0.2
z_max = 30.0
n_layers = int(z_max / layer_height)
print(f"モデル: 四角錐 底辺20mm 高さ30mm, レイヤー高さ {layer_height} mm")
print(f"総レイヤー数: {n_layers}")
print(f"{'z (mm)':>8} {'交線分数':>8} {'周長 (mm)':>12}")
total_path = 0.0
for z in [1.0, 6.0, 12.0, 18.0, 24.0, 29.0]:
    segs, perim = slice_at_z(pyr, z)
    print(f"{z:8.1f} {len(segs):8d} {perim:12.3f}")

# 全レイヤーの外周経路長を積算(印刷時間の概算に使用)
for i in range(1, n_layers + 1):
    z = i * layer_height
    _, perim = slice_at_z(pyr, min(z, z_max - 1e-6))
    total_path += perim
print_speed = 50.0  # mm/s
print(f"外周経路長の総和: {total_path:.1f} mm")
print(f"外周のみの印刷時間概算 ({print_speed:.0f} mm/s): {total_path/print_speed:.1f} s")

実行結果:

モデル: 四角錐 底辺20mm 高さ30mm, レイヤー高さ 0.2 mm
総レイヤー数: 150
  z (mm)     交線分数      周長 (mm)
     1.0        4       77.333
     6.0        4       64.000
    12.0        4       48.000
    18.0        4       32.000
    24.0        4       16.000
    29.0        4        2.667
外周経路長の総和: 5960.0 mm
外周のみの印刷時間概算 (50 mm/s): 119.2 s

各層で交線分がちょうど4本(正方形の4辺)得られ、周長が高さとともに線形に減少しています。z = 1 mmで周長77.3 mm、頂点近くのz = 29 mmでは2.67 mmと、四角錐の相似的な縮小を正しく捉えています。全150層の外周経路長を合計すると約5960 mmとなり、50 mm/sで外周のみを描くなら約119秒という概算が得られます。実際の造形ではこれにインフィルやトラベル移動、加減速が加わるため、これはあくまで下限の目安です。

💡 この実装の範囲と実務との差

ここでは「線分を集めて周長を測る」ところまでを実装しました。実用的なスライサーは、この線分群を端点でつないで閉じた輪郭(ループ)に整列させ、内外を判定してシェルとインフィルを生成し、さらにG-code(第1章)へ変換します。特に、複数の穴を持つ断面や自己交差する輪郭の処理、浮動小数点誤差で端点がわずかにずれた線分の接続は、実装上の難所です。本コードは幾何の核心を示すものであり、そのまま製品スライサーになるわけではありません。

5.3 熱シミュレーション:層の冷却

5.3.1 非定常熱伝導方程式

第2章では、層間接着が「界面がガラス転移温度(Tg)より高温である間だけ進む」ことを学びました。では、押し出された層は実際にどれくらいの速さで冷えるのでしょうか。これを定量的に扱うのが非定常熱伝導方程式(Transient Heat Conduction Equation)です。1次元では次の形になります。

∂T/∂t = α · ∂²T/∂x², α = k / (ρ·c)(熱拡散率)

ここで α は熱拡散率(thermal diffusivity)で、熱伝導率 k・密度 ρ・比熱 c から決まります。この偏微分方程式を、空間を格子に区切って差分で近似する陽解法(explicit finite difference)で解きます。陽解法は実装が明快ですが、時間刻み dt が大きすぎると数値的に発散するため、格子フーリエ数 r = α·dt/dx² が 0.5 以下という安定条件を守る必要があります。

ここでは、押出直後の厚さ0.8 mmの薄壁が、両面から周囲空気へ対流冷却されながら固化する過程を計算します。境界条件には対流(ニュートンの冷却則)をゴースト点法で組み込みます。

コード例3: 1次元非定常熱伝導の陽解法

import numpy as np

# 押出直後の薄壁が両面から対流冷却されながら固化する過程を
# 1次元非定常熱伝導方程式 dT/dt = alpha d2T/dx2 の陽解法で解く
L = 0.8e-3        # 壁厚 m (0.8 mm)
nx = 21           # 格子点数
dx = L / (nx - 1)
alpha = 1.3e-7    # PLAの熱拡散率 m^2/s
T_ext = 210.0     # 押出温度 C
T_env = 30.0      # 周囲温度 C
Tg = 60.0         # PLAのガラス転移温度 C
h = 40.0          # 対流熱伝達率 W/(m^2 K)
k = 0.13          # PLAの熱伝導率 W/(m K)

dt = 0.2 * dx**2 / alpha          # 安定条件 (Fourier数 <= 0.5) を満たす時間刻み
r = alpha * dt / dx**2
Bi = h * dx / k                   # 格子ビオ数
print(f"格子間隔 dx = {dx*1e6:.1f} um, 時間刻み dt = {dt*1000:.3f} ms")
print(f"格子フーリエ数 r = {r:.3f} (<=0.5 で安定), 格子ビオ数 Bi = {Bi:.4f}")

T = np.full(nx, T_ext)
t = 0.0
t_center_Tg = None
checkpoints = [0.0, 0.5, 1.0, 2.0, 5.0]
records = {}
next_cp = 0
max_steps = 2_000_000
for step in range(max_steps):
    Tn = T.copy()
    # 内部格子: 陽的差分
    Tn[1:-1] = T[1:-1] + r * (T[2:] - 2*T[1:-1] + T[:-2])
    # 対流境界(両端): ゴースト点法で対流を反映
    Tn[0] = T[0] + 2*r*(T[1] - T[0]) - 2*r*Bi*(T[0] - T_env)
    Tn[-1] = T[-1] + 2*r*(T[-2] - T[-1]) - 2*r*Bi*(T[-1] - T_env)
    T = Tn
    t += dt
    if t_center_Tg is None and T[nx//2] <= Tg:
        t_center_Tg = t
    while next_cp < len(checkpoints) and t >= checkpoints[next_cp]:
        records[checkpoints[next_cp]] = (T[0], T[nx//2])
        next_cp += 1
    if T.max() <= Tg and t_center_Tg is not None:
        break

print(f"{'時刻 (s)':>8} {'表面温度 (C)':>14} {'中心温度 (C)':>14}")
for cp in checkpoints:
    if cp in records:
        surf, cen = records[cp]
        print(f"{cp:8.1f} {surf:14.1f} {cen:14.1f}")
print(f"中心が Tg={Tg:.0f}C に達する時刻: {t_center_Tg:.2f} s")

実行結果:

格子間隔 dx = 40.0 um, 時間刻み dt = 2.462 ms
格子フーリエ数 r = 0.200 (<=0.5 で安定), 格子ビオ数 Bi = 0.0123
  時刻 (s)       表面温度 (C)       中心温度 (C)
     0.0          209.1          210.0
     0.5          194.7          204.8
     1.0          187.0          196.7
     2.0          172.6          181.5
     5.0          136.9          143.5
中心が Tg=60C に達する時刻: 18.86 s

格子フーリエ数 r = 0.20 は安定条件(≤0.5)を満たしており、計算は発散せず滑らかな冷却曲線を与えます。表面は空気と接するため中心よりわずかに速く冷えますが、格子ビオ数 Bi = 0.012 と非常に小さい(=内部の熱伝導が対流に対して十分速い)ため、表面と中心の温度差は小さく、壁全体がほぼ一様に冷えていきます。中心がTg(60°C)を下回るまで約19秒——これは周囲が30°Cと比較的暖かいためで、冷却ファンで周囲温度を下げれば、この時間は大きく短縮されます。第2章で学んだ「暖かい環境ほど溶着ウィンドウが長い」という定性的な結論が、ここでは具体的な秒数として裏づけられます。

⚠️ モデルの単純化

このモデルは1次元・対流冷却のみを扱い、下の層や造形プレートへの熱伝導、押出時の潜熱、温度依存の物性を無視しています。そのため実際の冷却より遅めに出る傾向があります。目的は絶対時間の精密予測ではなく、「安定条件を満たした陽解法で冷却曲線が得られること」「周囲温度と熱伝達が冷却の速さを支配すること」を体感することにあります。定量的な設計には、より高次の2次元・3次元モデルや有限要素法が用いられます。

5.4 プロセスパラメータの最適化

5.4.1 造形時間と表面品質のトレードオフ

3Dプリンティングの現場では、「できるだけ速く、できるだけきれいに」という相反する要求を両立させる必要があります。この2つはしばしば逆方向を向きます。レイヤー高さを大きくすれば層数が減って造形時間は短くなりますが、積層痕が目立って表面が粗くなります。ここでは、レイヤー高さと印刷速度という2つのパラメータをグリッドサーチ(grid search)で総当たりし、トレードオフを定量化します。

造形時間と表面粗さは、次の簡潔なモデルで見積もります。

そのうえで、「表面粗さ Ra が許容値6 µm以下」という品質制約を満たすなかで、造形時間が最短になる条件を探します。

コード例4: グリッドサーチによるパラメータ最適化

# 造形時間と表面品質のトレードオフをグリッドサーチで評価する。
# 造形時間 ~ 高さ/(レイヤー高さ) x 経路/(速度)
# 表面粗さ Ra ~ レイヤー高さに比例(積層痕)
part_height = 30.0      # mm
path_per_layer = 80.0   # mm/層 (外周+インフィルの代表値)
layer_heights = [0.10, 0.15, 0.20, 0.30]
speeds = [40, 60, 80, 100]     # mm/s
Ra_max = 6.0            # 許容表面粗さ um

print(f"{'LH(mm)':>7} {'v(mm/s)':>8} {'時間(min)':>10} {'Ra(um)':>8} {'許容':>6}")
results = []
for lh in layer_heights:
    n_layers = part_height / lh
    Ra = 1000.0 * lh**2 / 8.0    # 積層痕の幾何モデル Ra ~ LH^2/(8) を um換算
    for v in speeds:
        t_build = n_layers * path_per_layer / v / 60.0   # min
        ok = Ra <= Ra_max
        results.append((lh, v, t_build, Ra, ok))
        flag = "OK" if ok else "NG"
        print(f"{lh:7.2f} {v:8d} {t_build:10.2f} {Ra:8.2f} {flag:>6}")

feasible = [r for r in results if r[4]]
best = min(feasible, key=lambda x: x[2])
print(f"\n制約 Ra<={Ra_max:.0f}um を満たす最短時間の条件:")
print(f"  レイヤー高さ {best[0]:.2f} mm, 速度 {best[1]} mm/s, "
      f"造形時間 {best[2]:.2f} min, Ra {best[3]:.2f} um")

実行結果:

 LH(mm)  v(mm/s)    時間(min)   Ra(um)     許容
   0.10       40      10.00     1.25     OK
   0.10       60       6.67     1.25     OK
   0.10       80       5.00     1.25     OK
   0.10      100       4.00     1.25     OK
   0.15       40       6.67     2.81     OK
   0.15       60       4.44     2.81     OK
   0.15       80       3.33     2.81     OK
   0.15      100       2.67     2.81     OK
   0.20       40       5.00     5.00     OK
   0.20       60       3.33     5.00     OK
   0.20       80       2.50     5.00     OK
   0.20      100       2.00     5.00     OK
   0.30       40       3.33    11.25     NG
   0.30       60       2.22    11.25     NG
   0.30       80       1.67    11.25     NG
   0.30      100       1.33    11.25     NG

制約 Ra<=6um を満たす最短時間の条件:
  レイヤー高さ 0.20 mm, 速度 100 mm/s, 造形時間 2.00 min, Ra 5.00 um

レイヤー高さ0.30 mmはRa = 11.25 µmとなり品質制約を満たしません(NG)。造形時間だけを見れば0.30 mm・100 mm/sの1.33分が最速ですが、表面品質を犠牲にしています。品質制約を満たす範囲で最短なのはレイヤー高さ0.20 mm・速度100 mm/sで、造形時間2.00分・Ra 5.00 µmという結論が得られます。このように、複数の目的が競合するときは「制約付き最適化」として定式化すると、意思決定が明確になります。

💡 グリッドサーチの位置づけ

パラメータが2〜3個で、それぞれ数段階なら、総当たりのグリッドサーチが最も見通しがよく確実です。ただし、パラメータが増えると組み合わせ数が指数的に膨れ上がります(次元の呪い)。実務ではベイズ最適化や遺伝的アルゴリズムといった、より少ない試行で良い解を探す手法が使われます。ここで示したモデルの係数(Ra ≈ LH²/8 など)は説明のための単純化であり、実際の表面粗さはノズル形状・材料・冷却にも依存します。傾向をつかむ道具として活用してください。

5.5 機械学習によるプロセス予測

5.5.1 データ駆動でプロセスと品質を結ぶ

これまでのシミュレーションは、物理法則(幾何・熱伝導)を式で書き下す物理ベースモデルでした。一方、層間接着・欠陥・最終強度のように、多数の要因が複雑に絡んで単純な式にしづらい現象もあります。そこで有効なのが、実験・シミュレーションのデータからパラメータと結果の関係を学ぶ機械学習(Machine Learning)です。

ここでは、ノズル温度・レイヤー高さ・印刷速度・インフィル密度という4つのプロセスパラメータから引張強度を予測するモデルを作ります。実験データの代わりに、既知の生成関数(第2章の物理的傾向を反映:高温ほど溶着が進み強度増、厚い層ほど層間接着が低下、高速ほど溶着時間が減り強度低下、インフィルは強度にほぼ線形寄与)にノイズを加えた合成データ500点を用います。回帰器にはランダムフォレスト回帰(Random Forest Regression)を使い、決定係数 R²・平均絶対誤差(MAE)で精度を評価し、特徴量重要度から支配的なパラメータを読み取ります。

コード例5: ランダムフォレストによる強度予測

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_absolute_error

rng = np.random.default_rng(42)
n = 500
# 合成データ: プロセスパラメータ -> 引張強度(MPa)
nozzle_T = rng.uniform(190, 230, n)     # ノズル温度 C
layer_h  = rng.uniform(0.10, 0.30, n)   # レイヤー高さ mm
speed    = rng.uniform(40, 100, n)      # 印刷速度 mm/s
infill   = rng.uniform(20, 100, n)      # インフィル密度 %

# 既知の生成関数(層間溶着・インフィルの寄与)+ ノイズ
strength = (
    18.0
    + 0.22 * (nozzle_T - 190)          # 高温ほど層間溶着が進み強度増
    - 55.0 * (layer_h - 0.10)          # 厚い層ほど層間接着が低下
    - 0.05 * (speed - 40)              # 高速ほど溶着時間が減り強度低下
    + 0.28 * infill                    # インフィルは強度にほぼ線形寄与
    + rng.normal(0, 2.0, n)
)

X = np.column_stack([nozzle_T, layer_h, speed, infill])
y = strength
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=0)

model = RandomForestRegressor(n_estimators=200, random_state=0)
model.fit(X_tr, y_tr)
pred = model.predict(X_te)
r2 = r2_score(y_te, pred)
mae = mean_absolute_error(y_te, pred)

print(f"学習データ数: {len(X_tr)}, テストデータ数: {len(X_te)}")
print(f"テスト R^2 : {r2:.3f}")
print(f"テスト MAE : {mae:.2f} MPa")
names = ["ノズル温度", "レイヤー高さ", "印刷速度", "インフィル"]
print("特徴量重要度:")
for nm, imp in sorted(zip(names, model.feature_importances_), key=lambda x: -x[1]):
    print(f"  {nm:12s}: {imp:.3f}")

# 未知条件の強度予測
q = np.array([[220.0, 0.15, 50.0, 80.0]])
pq = model.predict(q)[0]
print(f"予測: ノズル220C/LH0.15mm/50mm/s/インフィル80% -> 引張強度 {pq:.1f} MPa")

実行結果:

学習データ数: 375, テストデータ数: 125
テスト R^2 : 0.906
テスト MAE : 1.93 MPa
特徴量重要度:
  インフィル       : 0.686
  レイヤー高さ      : 0.159
  ノズル温度       : 0.123
  印刷速度        : 0.031
予測: ノズル220C/LH0.15mm/50mm/s/インフィル80% -> 引張強度 42.7 MPa

テストデータでの決定係数 R² = 0.906、平均絶対誤差 1.93 MPaと、未知のパラメータ組合せに対しても良好な予測精度が得られました。特徴量重要度を見ると、インフィル密度(0.686)が圧倒的に支配的で、次いでレイヤー高さ・ノズル温度が続き、印刷速度の寄与は小さいことが分かります。これは生成関数でインフィルの係数(0.28×最大80)が最も大きく強度に効くように設定したことと整合しており、モデルが物理的傾向を正しく学習できたことを示します。実務では、こうした重要度分析が「どのパラメータを優先的に管理・最適化すべきか」の指針になります。

⚠️ 合成データであることの意味

ここでのデータは既知の関数から生成した合成データであり、モデルが高精度なのは「答えを知っている関数」を学習させたからです。実際の造形では、材料ロット差・装置個体差・環境揺らぎが加わり、同じR²は望めません。また、学習データの範囲外(外挿)では予測が信頼できない点にも注意が必要です。機械学習は万能の予測装置ではなく、十分な質と量のデータがあって初めて有効な道具です。ここでは「パイプラインの組み方」と「重要度の読み解き方」を学ぶことが目的です。

演習問題

理解を確認するための演習です。まず自分で考えてから解答を開いてください。

演習1(基礎): 水密性の意味

コード例1で、立方体メッシュの1枚の三角形をうっかり削除したとします。is_watertight の結果はどう変わり、それは造形にどんな問題を引き起こしますか。

解答を見る

三角形を1枚削除すると、その三角形が担っていた3本のエッジのうち、隣接三角形と共有していたエッジは共有数が2から1へ減ります。共有数が1のエッジが生じるため np.all(counts == 2) が False となり、watertight = False(穴あき)と判定されます。表面に穴が空いたメッシュは内外の区別が曖昧になり、スライサーが断面輪郭を閉じられず、造形が失敗したり意図しない中空・欠損が生じます。

演習2(計算): スライス周長の予測

コード例2の四角錐(底辺20 mm・高さ30 mm)を z = 15 mm でスライスすると、断面正方形の一辺と周長はいくらになりますか。相似比から手計算してください。

解答を見る

四角錐は高さ30 mmで頂点に収束するので、高さ z における断面正方形の一辺は 20 ×(1 − z/30) です。z = 15 mm では 20 ×(1 − 0.5) = 10 mm。周長は 4 × 10 = 40 mm。コード例2の出力(z=12で48、z=18で32)の中間として整合します。

演習3(考察): 安定条件

コード例3で、時間刻みを dt = 0.6 * dx**2 / alpha に変えると格子フーリエ数 r はいくつになり、計算結果はどうなると予想されますか。

解答を見る

r = α·dt/dx² = 0.6 となり、陽解法の安定条件 r ≤ 0.5 を超えます。このとき数値解は時間ステップごとに振動が増幅し、温度が物理的にありえない値(例えば数百度を超える発散や負の温度)へ暴走します。陽解法では安定条件が絶対で、dxを細かくするほどdtも小さくする必要があります。安定条件に縛られたくない場合は、無条件安定な陰解法(implicit method)を用います。

演習4(応用): 制約の変更

コード例4で許容表面粗さを Ra_max = 3.0(µm)に厳しくすると、最短造形時間の条件はどう変わりますか。表の値から答えてください。

解答を見る

Ra ≤ 3.0 µm を満たすのはレイヤー高さ0.10 mm(Ra=1.25)と0.15 mm(Ra=2.81)のみで、0.20 mm(Ra=5.00)は除外されます。この範囲で最短時間はレイヤー高さ0.15 mm・速度100 mm/sの2.67分です。品質要求を厳しくすると、選べるレイヤー高さの上限が下がり、造形時間が延びるというトレードオフが明確に現れます。

演習5(応用): 特徴量重要度の解釈

コード例5で、もし特徴量重要度において「印刷速度」が最大になったとしたら、生成関数や実データについてどんな解釈が可能ですか。

解答を見る

印刷速度の重要度が最大になるのは、強度が速度に強く依存する場合です。生成関数でいえば速度の係数(現状 −0.05)が他より大きい、あるいは速度が広い範囲で変動して強度への寄与が支配的になった状況です。実データでこれが起きたなら、その装置・材料では高速化に伴う溶着時間の不足(層間接着の低下)が強度のボトルネックになっている可能性を示唆します。対策として速度上限の見直しやノズル温度の引き上げ(第2章の層間強度モデル)が候補になります。重要度は「どの因子を優先して管理すべきか」を教えてくれる診断ツールです。

まとめ

本章では、シリーズで学んだ積層造形の知識を、5つの実行可能なPythonコードとして手を動かして確認しました。要点は次のとおりです。

✅ シリーズ完結

おめでとうございます。全5章を通じて、積層造形の原理(第1章)、材料押出法(第2章)、光造形と粉末床溶融結合(第3章)、材料噴射・結合剤噴射ほか(第4章)、そして本章のPython実践まで学び終えました。原理の理解と、それをコードで確かめる力の両方が身についたはずです。ここからは、実際の3Dプリンタやオープンソースのスライサー・メッシュライブラリ(trimesh など)に触れ、自分の課題に応用していってください。

次のステップ

本章で扱ったコードは、あえて外部ライブラリを最小限にして仕組みを露わにしたものです。実務や研究へ進むための発展として、次のような方向があります。

本シリーズが、あなたが積層造形を「使う」段階から「理解して設計・最適化する」段階へ進む足がかりとなれば幸いです。

参考文献

  1. Gibson, I., Rosen, D., & Stucker, B. (2021). Additive Manufacturing Technologies (3rd ed.). Springer. - AMのデータ処理・スライシング・シミュレーションを含む標準的教科書
  2. ISO/ASTM 52900:2021. Additive manufacturing — General principles — Fundamentals and vocabulary. - AMプロセス分類と用語の国際標準規格
  3. Incropera, F.P., DeWitt, D.P., Bergman, T.L., & Lavine, A.S. (2017). Fundamentals of Heat and Mass Transfer (8th ed.). Wiley. - 非定常熱伝導と陽解法・安定条件の標準的教科書
  4. Harris, C.R., et al. (2020). "Array Programming with NumPy." Nature, 585, 357-362. - 本章のメッシュ・数値計算を支えるNumPyの基盤論文
  5. Pedregosa, F., et al. (2011). "Scikit-learn: Machine Learning in Python." Journal of Machine Learning Research, 12, 2825-2830. - コード例5で用いたscikit-learnの基礎文献
  6. Dawson-Haggerty, M., et al. (2024). trimesh: Python library for loading and using triangular meshes. https://trimesh.org/ - 実STL処理・メッシュ修復・可視化のための包括的ライブラリ
  7. Slic3r Project. (2024). Slic3r Manual: Slicing Algorithms. https://manual.slic3r.org/ - 平面交差スライシングと輪郭生成の実装解説

使用ツールとライブラリ

免責事項