第4章：積層造形の基礎

flowchart LR
    A[1986
SLA発明
Chuck Hull] --> B[1988
SLS登場
Carl Deckard]
    B --> C[1992
FDM特許
Stratasys社]
    C --> D[2005
RepRap
オープンソース化]
    D --> E[2012
金属AM普及
EBM/SLM]
    E --> F[2023
産業化加速
大型・高速化]

    style A fill:#e3f2fd
    style B fill:#fff3e0
    style C fill:#e8f5e9
    style D fill:#f3e5f5
    style E fill:#fce4ec
    style F fill:#fff9c4
        

flowchart TD
    AM[積層造形
7つのプロセス] --> MEX[Material Extrusion
材料押出]
    AM --> VPP[Vat Photopolymerization
液槽光重合]
    AM --> PBF[Powder Bed Fusion
粉末床溶融結合]
    AM --> MJ[Material Jetting
材料噴射]
    AM --> BJ[Binder Jetting
結合剤噴射]
    AM --> SL[Sheet Lamination
シート積層]
    AM --> DED[Directed Energy Deposition
指向性エネルギー堆積]

    MEX --> MEX_EX[FDM/FFF
低コスト・普及型]
    VPP --> VPP_EX[SLA/DLP
高精度・高表面品質]
    PBF --> PBF_EX[SLS/SLM/EBM
高強度・金属対応]

    style AM fill:#f093fb
    style MEX fill:#e3f2fd
    style VPP fill:#fff3e0
    style PBF fill:#e8f5e9
    style MJ fill:#f3e5f5
    style BJ fill:#fce4ec
    style SL fill:#fff9c4
    style DED fill:#fce4ec
        

1.3 STLファイル形式とデータ処理

1.3.1 STLファイルの構造

STL（STereoLithography）は、AMで最も広く使用される3Dモデルファイル形式で、1987年に3D Systems社が開発しました。STLファイルは物体表面を三角形メッシュ（Triangle Mesh）の集合として表現します。

STLファイルの基本構造

ASCII STL形式の例：

solid cube
  facet normal 0 0 1
    outer loop
      vertex 0 0 10
      vertex 10 0 10
      vertex 10 10 10
    endloop
  endfacet
  facet normal 0 0 1
    outer loop
      vertex 0 0 10
      vertex 10 10 10
      vertex 0 10 10
    endloop
  endfacet
  ...
endsolid cube

STLフォーマットの2つの種類：

1. ASCII STL: 人間が読めるテキスト形式。ファイルサイズ大（同じモデルでBinaryの10-20倍）。デバッグ・検証に有用。
2. Binary STL: バイナリ形式、ファイルサイズ小、処理高速。産業用途で標準。構造：80バイトヘッダー + 4バイト（三角形数） + 各三角形50バイト（法線12B + 頂点36B + 属性2B）。

1.3.2 STLファイルの重要概念

1. 法線ベクトル（Normal Vector）

各三角形面には法線ベクトル（外向き方向）が定義され、物体の「内側」と「外側」を区別します。法線方向は右手の法則で決定されます：

頂点順序ルール： 頂点v1, v2, v3は反時計回り（CCW: Counter-ClockWise）に配置され、外から見て反時計回りの順序で法線が外向きになります。

2. 多様体（Manifold）条件

STLメッシュが3Dプリント可能であるためには、多様体（Manifold）でなければなりません：

• エッジ共有: すべてのエッジ（辺）は正確に2つの三角形に共有される
• 頂点共有: すべての頂点は連続した三角形扇（fan）に属する
• 閉じた表面: 穴や開口部がなく、完全に閉じた表面を形成
• 自己交差なし: 三角形が互いに交差・貫通していない

1.3.3 STLファイルの品質指標

STLメッシュの品質は以下の指標で評価されます：

1. 三角形数（Triangle Count）: 通常10,000-500,000個。過少（粗いモデル）または過多（ファイルサイズ大・処理遅延）は避ける。
2. エッジ長の一様性: 極端に大小の三角形が混在すると造形品質低下。理想的には0.1-1.0 mm範囲。
3. アスペクト比（Aspect Ratio）: 細長い三角形（高アスペクト比）は数値誤差の原因。理想的にはアスペクト比 < 10。
4. 法線の一貫性: すべての法線が外向き統一。反転法線が混在すると内外判定エラー。

1.3.4 Pythonライブラリによる STL処理

PythonでSTLファイルを扱うための主要ライブラリ：

1. numpy-stl: 高速STL読込・書込、体積・表面積計算、法線ベクトル操作。シンプルで軽量。
2. trimesh: 包括的な3Dメッシュ処理ライブラリ。メッシュ修復、ブーリアン演算、レイキャスト、衝突検出。多機能だが依存関係多い。
3. PyMesh: 高度なメッシュ処理（リメッシュ、サブディビジョン、フィーチャー抽出）。インストールやや複雑。

numpy-stlの基本的な使用法：

from stl import mesh
import numpy as np

# STLファイルを読み込み
your_mesh = mesh.Mesh.from_file('model.stl')

# 基本的な幾何情報
volume, cog, inertia = your_mesh.get_mass_properties()
print(f"Volume: {volume:.2f} mm³")
print(f"Center of Gravity: {cog}")
print(f"Surface Area: {your_mesh.areas.sum():.2f} mm²")

# 三角形数
print(f"Number of Triangles: {len(your_mesh.vectors)}")

1.4 スライシングとツールパス生成

STLファイルを3Dプリンタが理解できる指令（G-code）に変換するプロセスをスライシング（Slicing）といいます。このセクションでは、スライシングの基本原理、ツールパス戦略、そしてG-codeの基礎を学びます。

1.4.1 スライシングの基本原理

スライシングは、3Dモデルを一定の高さ（レイヤー高さ）で水平に切断し、各層の輪郭を抽出するプロセスです：

flowchart TD
    A[3Dモデル
STLファイル] --> B[Z軸方向に
層状にスライス]
    B --> C[各層の輪郭抽出
Contour Detection]
    C --> D[シェル生成
Perimeter Path]
    D --> E[インフィル生成
Infill Path]
    E --> F[サポート追加
Support Structure]
    F --> G[ツールパス最適化
Retraction/Travel]
    G --> H[G-code出力]

    style A fill:#e3f2fd
    style H fill:#e8f5e9
        

レイヤー高さ（Layer Height）の選択

レイヤー高さは造形品質と造形時間のトレードオフを決定する最重要パラメータです：

1.4.2 シェルとインフィル戦略

シェル（外殻）の生成

シェル（Shell/Perimeter）は、各層の外周部を形成する経路です：

• シェル数（Perimeter Count）: 通常2-4本。外部品質と強度に影響。
  • 1本: 非常に弱い、透明性高い、装飾用のみ
  • 2本: 標準（バランス良好）
  • 3-4本: 高強度、表面品質向上、気密性向上
• シェル順序: 内側→外側（Inside-Out）が一般的。外側→内側は表面品質重視時に使用。

インフィル（内部充填）パターン

インフィル（Infill）は内部構造を形成し、強度と材料使用量を制御します：

1.4.3 サポート構造の生成

オーバーハング角度が45度を超える部分は、サポート構造（Support Structure）が必要です：

サポートのタイプ

• Linear Support（直線サポート）: 垂直な柱状サポート。シンプルで除去しやすいが、材料使用量多い。
• Tree Support（ツリーサポート）: 樹木状に分岐するサポート。材料使用量30-50%削減、除去しやすい。CuraやPrusaSlicerで標準サポート。
• Interface Layers（接合層）: サポート上面に薄い接合層を設ける。除去しやすく、表面品質向上。通常2-4層。

サポート設定の重要パラメータ

1.4.4 G-codeの基礎

G-codeは、3DプリンタやCNCマシンを制御する標準的な数値制御言語です。各行が1つのコマンドを表します：

主要なG-codeコマンド

G-codeの例（造形開始部分）

; === Start G-code ===
M140 S60       ; ベッドを60°Cに加熱開始（非待機）
M104 S210      ; ノズルを210°Cに加熱開始（非待機）
G28            ; 全軸ホーミング
G29            ; オートレベリング（ベッドメッシュ計測）
M190 S60       ; ベッド温度到達を待機
M109 S210      ; ノズル温度到達を待機
G92 E0         ; 押出量をゼロリセット
G1 Z2.0 F3000  ; Z軸を2mm上昇（安全確保）
G1 X10 Y10 F5000  ; プライム位置へ移動
G1 Z0.3 F3000  ; Z軸を0.3mmへ降下（初層高さ）
G1 X100 E10 F1500 ; プライムライン描画（ノズル詰まり除去）
G92 E0         ; 押出量を再度ゼロリセット
; === 造形開始 ===

1.4.5 主要スライシングソフトウェア

1.4.6 ツールパス最適化戦略

効率的なツールパスは、造形時間・品質・材料使用量を改善します：

• リトラクション（Retraction）: 移動時にフィラメントを引き戻してストリング（糸引き）を防止。
  • 距離: 1-6mm（ボーデンチューブ式は4-6mm、ダイレクト式は1-2mm）
  • 速度: 25-45 mm/s
  • 過度なリトラクションはノズル詰まりの原因
• Z-hop（Z軸跳躍）: 移動時にノズルを上昇させて造形物との衝突を回避。0.2-0.5mm上昇。造形時間微増だが表面品質向上。
• コーミング（Combing）: 移動経路をインフィル上に制限し、表面への移動痕を低減。外観重視時に有効。
• シーム位置（Seam Position）: 各層の開始/終了点を揃える戦略。
  • Random: ランダム配置（目立たない）
  • Aligned: 一直線に配置（後加工でシームを除去しやすい）
  • Sharpest Corner: 最も鋭角なコーナーに配置（目立ちにくい）

Example 1: STLファイルの読み込みと基本情報取得

# ===================================
# Example 1: STLファイルの読み込みと基本情報取得
# ===================================

import numpy as np
from stl import mesh

# STLファイルを読み込む
your_mesh = mesh.Mesh.from_file('model.stl')

# 基本的な幾何情報を取得
volume, cog, inertia = your_mesh.get_mass_properties()

print("=== STLファイル基本情報 ===")
print(f"Volume: {volume:.2f} mm³")
print(f"Surface Area: {your_mesh.areas.sum():.2f} mm²")
print(f"Center of Gravity: [{cog[0]:.2f}, {cog[1]:.2f}, {cog[2]:.2f}] mm")
print(f"Number of Triangles: {len(your_mesh.vectors)}")

# バウンディングボックス（最小包含直方体）を計算
min_coords = your_mesh.vectors.min(axis=(0, 1))
max_coords = your_mesh.vectors.max(axis=(0, 1))
dimensions = max_coords - min_coords

print(f"\n=== バウンディングボックス ===")
print(f"X: {min_coords[0]:.2f} to {max_coords[0]:.2f} mm (幅: {dimensions[0]:.2f} mm)")
print(f"Y: {min_coords[1]:.2f} to {max_coords[1]:.2f} mm (奥行: {dimensions[1]:.2f} mm)")
print(f"Z: {min_coords[2]:.2f} to {max_coords[2]:.2f} mm (高さ: {dimensions[2]:.2f} mm)")

# 造形時間の簡易推定（レイヤー高さ0.2mm、速度50mm/sと仮定）
layer_height = 0.2  # mm
print_speed = 50    # mm/s
num_layers = int(dimensions[2] / layer_height)
# 簡易計算: 表面積に基づく推定
estimated_path_length = your_mesh.areas.sum() / layer_height  # mm
estimated_time_seconds = estimated_path_length / print_speed
estimated_time_minutes = estimated_time_seconds / 60

print(f"\n=== 造形推定 ===")
print(f"レイヤー数（0.2mm/層）: {num_layers} 層")
print(f"推定造形時間: {estimated_time_minutes:.1f} 分 ({estimated_time_minutes/60:.2f} 時間)")

# 出力例:
# === STLファイル基本情報 ===
# Volume: 12450.75 mm³
# Surface Area: 5832.42 mm²
# Center of Gravity: [25.34, 18.92, 15.67] mm
# Number of Triangles: 2456
#
# === バウンディングボックス ===
# X: 0.00 to 50.00 mm (幅: 50.00 mm)
# Y: 0.00 to 40.00 mm (奥行: 40.00 mm)
# Z: 0.00 to 30.00 mm (高さ: 30.00 mm)
#
# === 造形推定 ===
# レイヤー数（0.2mm/層）: 150 層
# 推定造形時間: 97.2 分 (1.62 時間)

Example 2: メッシュの法線ベクトル検証

# ===================================
# Example 2: メッシュの法線ベクトル検証
# ===================================

import numpy as np
from stl import mesh

def check_normals(mesh_data):
    """STLメッシュの法線ベクトルの整合性をチェック

    Args:
        mesh_data: numpy-stlのMeshオブジェクト

    Returns:
        tuple: (flipped_count, total_count, percentage)
    """
    # 右手系ルールで法線方向を確認
    flipped_count = 0
    total_count = len(mesh_data.vectors)

    for i, facet in enumerate(mesh_data.vectors):
        v0, v1, v2 = facet

        # エッジベクトルを計算
        edge1 = v1 - v0
        edge2 = v2 - v0

        # 外積で法線を計算（右手系）
        calculated_normal = np.cross(edge1, edge2)

        # 正規化
        norm = np.linalg.norm(calculated_normal)
        if norm > 1e-10:  # ゼロベクトルでないことを確認
            calculated_normal = calculated_normal / norm
        else:
            continue  # 縮退三角形をスキップ

        # ファイルに保存されている法線と比較
        stored_normal = mesh_data.normals[i]
        stored_norm = np.linalg.norm(stored_normal)

        if stored_norm > 1e-10:
            stored_normal = stored_normal / stored_norm

        # 内積で方向の一致をチェック
        dot_product = np.dot(calculated_normal, stored_normal)

        # 内積が負なら逆向き
        if dot_product < 0:
            flipped_count += 1

    percentage = (flipped_count / total_count) * 100 if total_count > 0 else 0

    return flipped_count, total_count, percentage

# STLファイルを読み込み
your_mesh = mesh.Mesh.from_file('model.stl')

# 法線チェックを実行
flipped, total, percent = check_normals(your_mesh)

print("=== 法線ベクトル検証結果 ===")
print(f"総三角形数: {total}")
print(f"反転法線数: {flipped}")
print(f"反転率: {percent:.2f}%")

if flipped == 0:
    print("\n✅ すべての法線が正しい方向を向いています")
    print("   このメッシュは3Dプリント可能です")
elif percent < 5:
    print("\n⚠️ 一部の法線が反転しています（軽微）")
    print("   スライサーが自動修正する可能性が高い")
else:
    print("\n❌ 多数の法線が反転しています（重大）")
    print("   メッシュ修復ツール（Meshmixer, netfabb）での修正を推奨")

# 出力例:
# === 法線ベクトル検証結果 ===
# 総三角形数: 2456
# 反転法線数: 0
# 反転率: 0.00%
#
# ✅ すべての法線が正しい方向を向いています
#    このメッシュは3Dプリント可能です

Example 3: マニフォールド性のチェック

# ===================================
# Example 3: マニフォールド性（Watertight）のチェック
# ===================================

import trimesh

# STLファイルを読み込み（trimeshは自動で修復を試みる）
mesh = trimesh.load('model.stl')

print("=== メッシュ品質診断 ===")

# 基本情報
print(f"Vertex count: {len(mesh.vertices)}")
print(f"Face count: {len(mesh.faces)}")
print(f"Volume: {mesh.volume:.2f} mm³")

# マニフォールド性をチェック
print(f"\n=== 3Dプリント適性チェック ===")
print(f"Is watertight (密閉性): {mesh.is_watertight}")
print(f"Is winding consistent (法線一致性): {mesh.is_winding_consistent}")
print(f"Is valid (幾何的妥当性): {mesh.is_valid}")

# 問題の詳細を診断
if not mesh.is_watertight:
    # 穴（hole）の数を検出
    try:
        edges = mesh.edges_unique
        edges_sorted = mesh.edges_sorted
        duplicate_edges = len(edges_sorted) - len(edges)
        print(f"\n⚠️ 問題検出:")
        print(f"   - メッシュに穴があります")
        print(f"   - 重複エッジ数: {duplicate_edges}")
    except:
        print(f"\n⚠️ メッシュ構造に問題があります")

# 修復を試みる
if not mesh.is_watertight or not mesh.is_winding_consistent:
    print(f"\n🔧 自動修復を実行中...")

    # 法線を修正
    trimesh.repair.fix_normals(mesh)
    print("   ✓ 法線ベクトルを修正")

    # 穴を埋める
    trimesh.repair.fill_holes(mesh)
    print("   ✓ 穴を充填")

    # 縮退三角形を削除
    mesh.remove_degenerate_faces()
    print("   ✓ 縮退面を削除")

    # 重複頂点を結合
    mesh.merge_vertices()
    print("   ✓ 重複頂点を結合")

    # 修復後の状態を確認
    print(f"\n=== 修復後の状態 ===")
    print(f"Is watertight: {mesh.is_watertight}")
    print(f"Is winding consistent: {mesh.is_winding_consistent}")

    # 修復したメッシュを保存
    if mesh.is_watertight:
        mesh.export('model_repaired.stl')
        print(f"\n✅ 修復完了！ model_repaired.stl として保存しました")
    else:
        print(f"\n❌ 自動修復失敗。Meshmixer等の専用ツールを推奨")
else:
    print(f"\n✅ このメッシュは3Dプリント可能です")

# 出力例:
# === メッシュ品質診断 ===
# Vertex count: 1534
# Face count: 2456
# Volume: 12450.75 mm³
#
# === 3Dプリント適性チェック ===
# Is watertight (密閉性): True
# Is winding consistent (法線一致性): True
# Is valid (幾何的妥当性): True
#
# ✅ このメッシュは3Dプリント可能です

Example 4: 基本的なスライシングアルゴリズム

# ===================================
# Example 4: 基本的なスライシングアルゴリズム
# ===================================

import numpy as np
from stl import mesh

def slice_mesh_at_height(mesh_data, z_height):
    """温度プロファイルを生成

    Args:
        t (array): 時間配列 [min]
        T_target (float): 保持温度 [°C]
        heating_rate (float): 加熱速度 [°C/min]
        hold_time (float): 保持時間 [min]
        cooling_rate (float): 冷却速度 [°C/min]

    Returns:
        array: 温度プロファイル [°C]
    """
    T_room = 25  # 室温
    T = np.zeros_like(t)

    # 加熱時間
    t_heat = (T_target - T_room) / heating_rate

    # 冷却開始時刻
    t_cool_start = t_heat + hold_time

    for i, time in enumerate(t):
        if time <= t_heat:
            # 加熱フェーズ
            T[i] = T_room + heating_rate * time
        elif time <= t_cool_start:
            # 保持フェーズ
            T[i] = T_target
        else:
            # 冷却フェーズ
            T[i] = T_target - cooling_rate * (time - t_cool_start)
            T[i] = max(T[i], T_room)  # 室温以下にはならない

    return T

def simulate_reaction_progress(T, t, Ea, D0, r0):
    """温度プロファイルに基づく反応進行を計算

    Args:
        T (array): 温度プロファイル [°C]
        t (array): 時間配列 [min]
        Ea (float): 活性化エネルギー [J/mol]
        D0 (float): 頻度因子 [m²/s]
        r0 (float): 粒子半径 [m]

    Returns:
        array: 反応率
    """
    R = 8.314
    C0 = 10000
    alpha = np.zeros_like(t)

    for i in range(1, len(t)):
        T_k = T[i] + 273.15
        D = D0 * np.exp(-Ea / (R * T_k))
        k = D * C0 / r0**2

        dt = (t[i] - t[i-1]) * 60  # min → s

        # 簡易積分（微小時間での反応進行）
        if alpha[i-1] < 0.99:
            dalpha = k * dt / (2 * (1 - (1-alpha[i-1])**(1/3)))
            alpha[i] = min(alpha[i-1] + dalpha, 1.0)
        else:
            alpha[i] = alpha[i-1]

    return alpha

# パラメータ設定
T_target = 1200  # °C
hold_time = 240  # min (4 hours)
Ea = 300e3  # J/mol
D0 = 5e-4  # m²/s
r0 = 5e-6  # m

# 異なる加熱速度での比較
heating_rates = [2, 5, 10, 20]  # °C/min
cooling_rate = 3  # °C/min

# 時間配列
t_max = 800  # min
t = np.linspace(0, t_max, 2000)

# プロット
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))

# 温度プロファイル
for hr in heating_rates:
    T_profile = temperature_profile(t, T_target, hr, hold_time, cooling_rate)
    ax1.plot(t/60, T_profile, linewidth=2, label=f'{hr}°C/min')

ax1.set_xlabel('Time (hours)', fontsize=12)
ax1.set_ylabel('Temperature (°C)', fontsize=12)
ax1.set_title('Temperature Profiles', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0, t_max/60])

# 反応進行
for hr in heating_rates:
    T_profile = temperature_profile(t, T_target, hr, hold_time, cooling_rate)
    alpha = simulate_reaction_progress(T_profile, t, Ea, D0, r0)
    ax2.plot(t/60, alpha, linewidth=2, label=f'{hr}°C/min')

ax2.axhline(y=0.95, color='red', linestyle='--', linewidth=1, label='Target (95%)')
ax2.set_xlabel('Time (hours)', fontsize=12)
ax2.set_ylabel('Conversion', fontsize=12)
ax2.set_title('Reaction Progress', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim([0, t_max/60])
ax2.set_ylim([0, 1])

plt.tight_layout()
plt.savefig('temperature_profile_optimization.png', dpi=300, bbox_inches='tight')
plt.show()

# 各加熱速度での95%反応到達時間を計算
print("\n95%反応到達時間の比較:")
print("=" * 60)
for hr in heating_rates:
    T_profile = temperature_profile(t, T_target, hr, hold_time, cooling_rate)
    alpha = simulate_reaction_progress(T_profile, t, Ea, D0, r0)

    # 95%到達時刻
    idx_95 = np.where(alpha >= 0.95)[0]
    if len(idx_95) > 0:
        t_95 = t[idx_95[0]] / 60
        print(f"加熱速度 {hr:2d}°C/min: t₉₅ = {t_95:.1f} hours")
    else:
        print(f"加熱速度 {hr:2d}°C/min: 反応不完全")

# 出力例:
# 95%反応到達時間の比較:
# ============================================================
# 加熱速度  2°C/min: t₉₅ = 7.8 hours
# 加熱速度  5°C/min: t₉₅ = 7.2 hours
# 加熱速度 10°C/min: t₉₅ = 6.9 hours
# 加熱速度 20°C/min: t₉₅ = 6.7 hours

演習問題

1.5.1 pycalphadとは

pycalphadは、CALPHAD（CALculation of PHAse Diagrams）法に基づく相図計算のためのPythonライブラリです。熱力学データベースから平衡相を計算し、反応経路の設計に有用です。

1.5.2 二元系相図の計算例

# ===================================
# Example 5: pycalphadで相図計算
# ===================================

# 注意: pycalphadのインストールが必要
# pip install pycalphad

from pycalphad import Database, equilibrium, variables as v
import matplotlib.pyplot as plt
import numpy as np

# TDBデータベースを読み込み（ここでは簡易的な例）
# 実際には適切なTDBファイルが必要
# 例: BaO-TiO2系

# 簡易的なTDB文字列（実際はより複雑）
tdb_string = """
$ BaO-TiO2 system (simplified)
ELEMENT BA   BCC_A2  137.327   !
ELEMENT TI   HCP_A3   47.867   !
ELEMENT O    GAS      15.999   !

FUNCTION GBCCBA   298.15  +GHSERBA;   6000 N !
FUNCTION GHCPTI   298.15  +GHSERTI;   6000 N !
FUNCTION GGASO    298.15  +GHSERO;    6000 N !

PHASE LIQUID:L %  1  1.0  !
PHASE BAO_CUBIC %  2  1 1  !
PHASE TIO2_RUTILE %  2  1 2  !
PHASE BATIO3 %  3  1 1 3  !
"""

# 注: 実際の計算には正式なTDBファイルが必要
# ここでは概念的な説明に留める

print("pycalphadによる相図計算の概念:")
print("=" * 60)
print("1. TDBデータベース（熱力学データ）を読み込む")
print("2. 温度・組成範囲を設定")
print("3. 平衡計算を実行")
print("4. 安定相を可視化")
print()
print("実際の適用例:")
print("- BaO-TiO2系: BaTiO3の形成温度・組成範囲")
print("- Si-N系: Si3N4の安定領域")
print("- 多元系セラミックスの相関係")

# 概念的なプロット（実データに基づくイメージ）
fig, ax = plt.subplots(figsize=(10, 7))

# 温度範囲
T = np.linspace(800, 1600, 100)

# 各相の安定領域（概念図）
# BaO + TiO2 → BaTiO3 反応
BaO_region = np.ones_like(T) * 0.3
TiO2_region = np.ones_like(T) * 0.7
BaTiO3_region = np.where((T > 1100) & (T < 1400), 0.5, np.nan)

ax.fill_between(T, 0, BaO_region, alpha=0.3, color='blue', label='BaO + TiO2')
ax.fill_between(T, BaO_region, TiO2_region, alpha=0.3, color='green',
                label='BaTiO3 stable')
ax.fill_between(T, TiO2_region, 1, alpha=0.3, color='red', label='Liquid')

ax.axhline(y=0.5, color='black', linestyle='--', linewidth=2,
           label='BaTiO3 composition')
ax.axvline(x=1100, color='gray', linestyle=':', linewidth=1, alpha=0.5)
ax.axvline(x=1400, color='gray', linestyle=':', linewidth=1, alpha=0.5)

ax.set_xlabel('Temperature (°C)', fontsize=12)
ax.set_ylabel('Composition (BaO mole fraction)', fontsize=12)
ax.set_title('Conceptual Phase Diagram: BaO-TiO2', fontsize=14, fontweight='bold')
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_xlim([800, 1600])
ax.set_ylim([0, 1])

# テキスト注釈
ax.text(1250, 0.5, 'BaTiO₃\nformation\nregion',
        fontsize=11, ha='center', va='center',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.7))

plt.tight_layout()
plt.savefig('phase_diagram_concept.png', dpi=300, bbox_inches='tight')
plt.show()

# 実際の使用例（コメントアウト）
"""
# 実際のpycalphad使用例
db = Database('BaO-TiO2.tdb')  # TDBファイル読み込み

# 平衡計算
eq = equilibrium(db, ['BA', 'TI', 'O'], ['LIQUID', 'BATIO3'],
                 {v.X('BA'): (0, 1, 0.01),
                  v.T: (1000, 1600, 50),
                  v.P: 101325})

# 結果プロット
eq.plot()
"""

1.6 実験計画法（DOE）による条件最適化

1.6.1 DOEとは

実験計画法（Design of Experiments, DOE）は、複数のパラメータが相互作用する系で、最小の実験回数で最適条件を見つける統計手法です。

固相反応で最適化すべき主要パラメータ：

• 反応温度（T）
• 保持時間（t）
• 粒子サイズ（r）
• 原料比（モル比）
• 雰囲気（空気、窒素、真空など）

1.6.2 応答曲面法（Response Surface Methodology）

# ===================================
# Example 6: DOEによる条件最適化
# ===================================

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize

# 仮想的な反応率モデル（温度と時間の関数）
def reaction_yield(T, t, noise=0):
    """温度と時間から反応率を計算（仮想モデル）

    Args:
        T (float): 温度 [°C]
        t (float): 時間 [hours]
        noise (float): ノイズレベル

    Returns:
        float: 反応率 [%]
    """
    # 最適値: T=1200°C, t=6 hours
    T_opt = 1200
    t_opt = 6

    # 二次モデル（ガウス型）
    yield_val = 100 * np.exp(-((T-T_opt)/150)**2 - ((t-t_opt)/3)**2)

    # ノイズ追加
    if noise > 0:
        yield_val += np.random.normal(0, noise)

    return np.clip(yield_val, 0, 100)

# 実験点配置（中心複合計画法）
T_levels = [1000, 1100, 1200, 1300, 1400]  # °C
t_levels = [2, 4, 6, 8, 10]  # hours

# グリッドで実験点を配置
T_grid, t_grid = np.meshgrid(T_levels, t_levels)
yield_grid = np.zeros_like(T_grid, dtype=float)

# 各実験点で反応率を測定（シミュレーション）
for i in range(len(t_levels)):
    for j in range(len(T_levels)):
        yield_grid[i, j] = reaction_yield(T_grid[i, j], t_grid[i, j], noise=2)

# 結果の表示
print("実験計画法による反応条件最適化")
print("=" * 70)
print(f"{'Temperature (°C)':<20} {'Time (hours)':<20} {'Yield (%)':<20}")
print("-" * 70)
for i in range(len(t_levels)):
    for j in range(len(T_levels)):
        print(f"{T_grid[i, j]:<20} {t_grid[i, j]:<20} {yield_grid[i, j]:<20.1f}")

# 最大反応率の条件を探す
max_idx = np.unravel_index(np.argmax(yield_grid), yield_grid.shape)
T_best = T_grid[max_idx]
t_best = t_grid[max_idx]
yield_best = yield_grid[max_idx]

print("-" * 70)
print(f"最適条件: T = {T_best}°C, t = {t_best} hours")
print(f"最大反応率: {yield_best:.1f}%")

# 3Dプロット
fig = plt.figure(figsize=(14, 6))

# 3D表面プロット
ax1 = fig.add_subplot(121, projection='3d')
T_fine = np.linspace(1000, 1400, 50)
t_fine = np.linspace(2, 10, 50)
T_mesh, t_mesh = np.meshgrid(T_fine, t_fine)
yield_mesh = np.zeros_like(T_mesh)

for i in range(len(t_fine)):
    for j in range(len(T_fine)):
        yield_mesh[i, j] = reaction_yield(T_mesh[i, j], t_mesh[i, j])

surf = ax1.plot_surface(T_mesh, t_mesh, yield_mesh, cmap='viridis',
                        alpha=0.8, edgecolor='none')
ax1.scatter(T_grid, t_grid, yield_grid, color='red', s=50,
            label='Experimental points')

ax1.set_xlabel('Temperature (°C)', fontsize=10)
ax1.set_ylabel('Time (hours)', fontsize=10)
ax1.set_zlabel('Yield (%)', fontsize=10)
ax1.set_title('Response Surface', fontsize=12, fontweight='bold')
ax1.view_init(elev=25, azim=45)
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=5)

# 等高線プロット
ax2 = fig.add_subplot(122)
contour = ax2.contourf(T_mesh, t_mesh, yield_mesh, levels=20, cmap='viridis')
ax2.contour(T_mesh, t_mesh, yield_mesh, levels=10, colors='black',
            alpha=0.3, linewidths=0.5)
ax2.scatter(T_grid, t_grid, c=yield_grid, s=100, edgecolors='red',
            linewidths=2, cmap='viridis')
ax2.scatter(T_best, t_best, color='red', s=300, marker='*',
            edgecolors='white', linewidths=2, label='Optimum')

ax2.set_xlabel('Temperature (°C)', fontsize=11)
ax2.set_ylabel('Time (hours)', fontsize=11)
ax2.set_title('Contour Map', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
fig.colorbar(contour, ax=ax2, label='Yield (%)')

plt.tight_layout()
plt.savefig('doe_optimization.png', dpi=300, bbox_inches='tight')
plt.show()

1.6.3 実験計画の実践的アプローチ

実際の固相反応では、以下の手順でDOEを適用します：

1. スクリーニング実験（2水準要因計画法）: 影響の大きいパラメータを特定
2. 応答曲面法（中心複合計画法）: 最適条件の探索
3. 確認実験: 予測された最適条件で実験し、モデルを検証

1.7 反応速度曲線のフィッティング

1.7.1 実験データからの速度定数決定

# ===================================
# Example 7: 反応速度曲線フィッティング
# ===================================

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 実験データ（時間 vs 反応率）
# 例: BaTiO3合成 @ 1200°C
time_exp = np.array([0, 1, 2, 3, 4, 6, 8, 10, 12, 15, 20])  # hours
conversion_exp = np.array([0, 0.15, 0.28, 0.38, 0.47, 0.60,
                          0.70, 0.78, 0.84, 0.90, 0.95])

# Jander式モデル
def jander_model(t, k):
    """Jander式による反応率計算

    Args:
        t (array): 時間 [hours]
        k (float): 速度定数

    Returns:
        array: 反応率
    """
    # [1 - (1-α)^(1/3)]² = kt を α について解く
    kt = k * t
    alpha = 1 - (1 - np.sqrt(kt))**3
    alpha = np.clip(alpha, 0, 1)  # 0-1の範囲に制限
    return alpha

# Ginstling-Brounshtein式（別の拡散モデル）
def gb_model(t, k):
    """Ginstling-Brounshtein式

    Args:
        t (array): 時間
        k (float): 速度定数

    Returns:
        array: 反応率
    """
    # 1 - 2α/3 - (1-α)^(2/3) = kt
    # 数値的に解く必要があるが、ここでは近似式を使用
    kt = k * t
    alpha = 1 - (1 - kt/2)**(3/2)
    alpha = np.clip(alpha, 0, 1)
    return alpha

# Power law (経験式)
def power_law_model(t, k, n):
    """べき乗則モデル

    Args:
        t (array): 時間
        k (float): 速度定数
        n (float): 指数

    Returns:
        array: 反応率
    """
    alpha = k * t**n
    alpha = np.clip(alpha, 0, 1)
    return alpha

# 各モデルでフィッティング
# Jander式
popt_jander, _ = curve_fit(jander_model, time_exp, conversion_exp, p0=[0.01])
k_jander = popt_jander[0]

# Ginstling-Brounshtein式
popt_gb, _ = curve_fit(gb_model, time_exp, conversion_exp, p0=[0.01])
k_gb = popt_gb[0]

# Power law
popt_power, _ = curve_fit(power_law_model, time_exp, conversion_exp, p0=[0.1, 0.5])
k_power, n_power = popt_power

# 予測曲線生成
t_fit = np.linspace(0, 20, 200)
alpha_jander = jander_model(t_fit, k_jander)
alpha_gb = gb_model(t_fit, k_gb)
alpha_power = power_law_model(t_fit, k_power, n_power)

# 残差計算
residuals_jander = conversion_exp - jander_model(time_exp, k_jander)
residuals_gb = conversion_exp - gb_model(time_exp, k_gb)
residuals_power = conversion_exp - power_law_model(time_exp, k_power, n_power)

# R²計算
def r_squared(y_true, y_pred):
    ss_res = np.sum((y_true - y_pred)**2)
    ss_tot = np.sum((y_true - np.mean(y_true))**2)
    return 1 - (ss_res / ss_tot)

r2_jander = r_squared(conversion_exp, jander_model(time_exp, k_jander))
r2_gb = r_squared(conversion_exp, gb_model(time_exp, k_gb))
r2_power = r_squared(conversion_exp, power_law_model(time_exp, k_power, n_power))

# プロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# フィッティング結果
ax1.plot(time_exp, conversion_exp, 'ko', markersize=8, label='Experimental data')
ax1.plot(t_fit, alpha_jander, 'b-', linewidth=2,
         label=f'Jander (R²={r2_jander:.4f})')
ax1.plot(t_fit, alpha_gb, 'r-', linewidth=2,
         label=f'Ginstling-Brounshtein (R²={r2_gb:.4f})')
ax1.plot(t_fit, alpha_power, 'g-', linewidth=2,
         label=f'Power law (R²={r2_power:.4f})')

ax1.set_xlabel('Time (hours)', fontsize=12)
ax1.set_ylabel('Conversion', fontsize=12)
ax1.set_title('Kinetic Model Fitting', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0, 20])
ax1.set_ylim([0, 1])

# 残差プロット
ax2.plot(time_exp, residuals_jander, 'bo-', label='Jander')
ax2.plot(time_exp, residuals_gb, 'ro-', label='Ginstling-Brounshtein')
ax2.plot(time_exp, residuals_power, 'go-', label='Power law')
ax2.axhline(y=0, color='black', linestyle='--', linewidth=1)

ax2.set_xlabel('Time (hours)', fontsize=12)
ax2.set_ylabel('Residuals', fontsize=12)
ax2.set_title('Residual Plot', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('kinetic_fitting.png', dpi=300, bbox_inches='tight')
plt.show()

# 結果サマリー
print("\n反応速度モデルのフィッティング結果:")
print("=" * 70)
print(f"{'Model':<25} {'Parameter':<30} {'R²':<10}")
print("-" * 70)
print(f"{'Jander':<25} {'k = ' + f'{k_jander:.4f} h⁻¹':<30} {r2_jander:.4f}")
print(f"{'Ginstling-Brounshtein':<25} {'k = ' + f'{k_gb:.4f} h⁻¹':<30} {r2_gb:.4f}")
print(f"{'Power law':<25} {'k = ' + f'{k_power:.4f}, n = {n_power:.4f}':<30} {r2_power:.4f}")
print("=" * 70)
print(f"\n最適モデル: {'Jander' if r2_jander == max(r2_jander, r2_gb, r2_power) else 'GB' if r2_gb == max(r2_jander, r2_gb, r2_power) else 'Power law'}")

# 出力例:
# 反応速度モデルのフィッティング結果:
# ======================================================================
# Model                     Parameter                      R²
# ----------------------------------------------------------------------
# Jander                    k = 0.0289 h⁻¹                 0.9953
# Ginstling-Brounshtein     k = 0.0412 h⁻¹                 0.9867
# Power law                 k = 0.2156, n = 0.5234         0.9982
# ======================================================================
#
# 最適モデル: Power law

1.8 高度なトピック: 微細構造制御

1.8.1 粒成長の抑制

固相反応では、高温・長時間保持により望ましくない粒成長が起こります。これを抑制する戦略：

• Two-step sintering: 高温で短時間保持後、低温で長時間保持
• 添加剤の使用: 粒成長抑制剤（例: MgO, Al₂O₃）を微量添加
• Spark Plasma Sintering (SPS): 急速加熱・短時間焼結

1.8.2 反応の機械化学的活性化

メカノケミカル法（高エネルギーボールミル）により、固相反応を室温付近で進行させることも可能です：

# ===================================
# Example 8: 粒成長シミュレーション
# ===================================

import numpy as np
import matplotlib.pyplot as plt

def grain_growth(t, T, D0, Ea, G0, n):
    """粒成長の時間発展

    Burke-Turnbull式: G^n - G0^n = k*t

    Args:
        t (array): 時間 [hours]
        T (float): 温度 [K]
        D0 (float): 頻度因子
        Ea (float): 活性化エネルギー [J/mol]
        G0 (float): 初期粒径 [μm]
        n (float): 粒成長指数（通常2-4）

    Returns:
        array: 粒径 [μm]
    """
    R = 8.314
    k = D0 * np.exp(-Ea / (R * T))
    G = (G0**n + k * t * 3600)**(1/n)  # hours → seconds
    return G

# パラメータ設定
D0_grain = 1e8  # μm^n/s
Ea_grain = 400e3  # J/mol
G0 = 0.5  # μm
n = 3

# 温度の影響
temps_celsius = [1100, 1200, 1300]
t_range = np.linspace(0, 12, 100)  # 0-12 hours

plt.figure(figsize=(12, 5))

# 温度依存性
plt.subplot(1, 2, 1)
for T_c in temps_celsius:
    T_k = T_c + 273.15
    G = grain_growth(t_range, T_k, D0_grain, Ea_grain, G0, n)
    plt.plot(t_range, G, linewidth=2, label=f'{T_c}°C')

plt.axhline(y=1.0, color='red', linestyle='--', linewidth=1,
            label='Target grain size')
plt.xlabel('Time (hours)', fontsize=12)
plt.ylabel('Grain Size (μm)', fontsize=12)
plt.title('Grain Growth at Different Temperatures', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.ylim([0, 5])

# Two-step sinteringの効果
plt.subplot(1, 2, 2)

# Conventional sintering: 1300°C, 6 hours
t_conv = np.linspace(0, 6, 100)
T_conv = 1300 + 273.15
G_conv = grain_growth(t_conv, T_conv, D0_grain, Ea_grain, G0, n)

# Two-step: 1300°C 1h → 1200°C 5h
t1 = np.linspace(0, 1, 20)
G1 = grain_growth(t1, 1300+273.15, D0_grain, Ea_grain, G0, n)
G_intermediate = G1[-1]

t2 = np.linspace(0, 5, 80)
G2 = grain_growth(t2, 1200+273.15, D0_grain, Ea_grain, G_intermediate, n)

t_two_step = np.concatenate([t1, t2 + 1])
G_two_step = np.concatenate([G1, G2])

plt.plot(t_conv, G_conv, 'r-', linewidth=2, label='Conventional (1300°C)')
plt.plot(t_two_step, G_two_step, 'b-', linewidth=2, label='Two-step (1300°C→1200°C)')
plt.axvline(x=1, color='gray', linestyle=':', linewidth=1, alpha=0.5)

plt.xlabel('Time (hours)', fontsize=12)
plt.ylabel('Grain Size (μm)', fontsize=12)
plt.title('Two-Step Sintering Strategy', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.ylim([0, 5])

plt.tight_layout()
plt.savefig('grain_growth_control.png', dpi=300, bbox_inches='tight')
plt.show()

# 最終粒径の比較
G_final_conv = grain_growth(6, 1300+273.15, D0_grain, Ea_grain, G0, n)
G_final_two_step = G_two_step[-1]

print("\n粒成長の比較:")
print("=" * 50)
print(f"Conventional (1300°C, 6h): {G_final_conv:.2f} μm")
print(f"Two-step (1300°C 1h + 1200°C 5h): {G_final_two_step:.2f} μm")
print(f"粒径抑制効果: {(1 - G_final_two_step/G_final_conv)*100:.1f}%")

# 出力例:
# 粒成長の比較:
# ==================================================
# Conventional (1300°C, 6h): 4.23 μm
# Two-step (1300°C 1h + 1200°C 5h): 2.87 μm
# 粒径抑制効果: 32.2%

学習目標の確認

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

基本理解

• ✅ 固相反応の3つの律速段階（核生成・界面反応・拡散）を説明できる
• ✅ Arrhenius式の物理的意味と温度依存性を理解している
• ✅ Jander式とGinstling-Brounshtein式の違いを説明できる
• ✅ 温度プロファイルの3要素（加熱速度・保持時間・冷却速度）の重要性を理解している

実践スキル

• ✅ Pythonで拡散係数の温度依存性をシミュレートできる
• ✅ Jander式を用いて反応進行を予測できる
• ✅ Kissinger法でDSC/TGデータから活性化エネルギーを計算できる
• ✅ DOE（実験計画法）で反応条件を最適化できる
• ✅ pycalphadを用いた相図計算の基礎を理解している

応用力

• ✅ 新規セラミックス材料の合成プロセスを設計できる
• ✅ 実験データから反応機構を推定し、適切な速度式を選択できる
• ✅ 産業プロセスでの条件最適化戦略を立案できる
• ✅ 粒成長制御の戦略（Two-step sintering等）を提案できる

演習問題

Easy（基礎確認）

Medium（応用）

import trimesh

mesh = trimesh.load('model.stl')

# ここにコードを追加: マニフォールド性をチェックし、
# 問題があれば自動修復を行い、修復後のメッシュを
# 'model_fixed.stl'として保存してください

import trimesh

mesh = trimesh.load('model.stl')

# マニフォールド性をチェック
print(f"Is watertight: {mesh.is_watertight}")
print(f"Is winding consistent: {mesh.is_winding_consistent}")

# 問題がある場合は修復
if not mesh.is_watertight or not mesh.is_winding_consistent:
    print("メッシュ修復を実行中...")

    # 法線を修正
    trimesh.repair.fix_normals(mesh)

    # 穴を埋める
    trimesh.repair.fill_holes(mesh)

    # 縮退三角形を削除
    mesh.remove_degenerate_faces()

    # 重複頂点を結合
    mesh.merge_vertices()

    # 修復結果を確認
    print(f"修復後 watertight: {mesh.is_watertight}")

    # 修復したメッシュを保存
    if mesh.is_watertight:
        mesh.export('model_fixed.stl')
        print("修復完了: model_fixed.stl として保存")
    else:
        print("⚠️ 自動修復失敗。Meshmixer等を使用してください")
else:
    print("✓ メッシュは3Dプリント可能です")

DOEの利点（従来法との比較）:

1. 交互作用の検出が可能
   • 従来法: 温度の影響、時間の影響を個別に評価
   • DOE: 「高温では時間を短くできる」といった交互作用を定量化
   • 例: 1300°Cでは4時間で十分だが、1100°Cでは8時間必要、など
2. 実験回数の削減
   • 従来法（OFAT: One Factor At a Time）:
     • 温度検討: 3回（時間固定）
     • 時間検討: 3回（温度固定）
     • 確認実験: 複数回
     • 合計: 10回以上
   • DOE: 9回で完了（全条件網羅＋交互作用解析）
   • さらに中心複合計画法を使えば7回に削減可能

追加の利点:

• 統計的に有意な結論が得られる（誤差評価が可能）
• 応答曲面を構築でき、未実施条件の予測が可能
• 最適条件が実験範囲外にある場合でも検出できる

Hard（発展）

次のステップ

第4章では積層造形（AM）の基礎として、ISO/ASTM 52900による7つのプロセス分類、STLファイル形式の構造、スライシングとG-codeの基本を学びました。次の第2章では、材料押出（FDM/FFF）の詳細な造形プロセス、材料特性、プロセスパラメータ最適化について学びます。

参考文献

1. Gibson, I., Rosen, D., & Stucker, B. (2015). Additive Manufacturing Technologies: 3D Printing, Rapid Prototyping, and Direct Digital Manufacturing (2nd ed.). Springer. pp. 1-35, 89-145, 287-334. - AM技術の包括的教科書、7つのプロセスカテゴリとSTLデータ処理の詳細解説
2. ISO/ASTM 52900:2021. Additive manufacturing — General principles — Fundamentals and vocabulary. International Organization for Standardization. - AM用語とプロセス分類の国際標準規格、産業界で広く参照される
3. Kruth, J.P., Leu, M.C., & Nakagawa, T. (1998). "Progress in Additive Manufacturing and Rapid Prototyping." CIRP Annals - Manufacturing Technology, 47(2), 525-540. - 選択的レーザー焼結とバインディング機構の理論的基礎
4. Hull, C.W. (1986). Apparatus for production of three-dimensional objects by stereolithography. US Patent 4,575,330. - 世界初のAM技術（SLA）の特許、AM産業の起源となる重要文献
5. Wohlers, T. (2023). Wohlers Report 2023: 3D Printing and Additive Manufacturing Global State of the Industry. Wohlers Associates, Inc. pp. 15-89, 156-234. - AM市場動向と産業応用の最新統計レポート、年次更新される業界標準資料
6. 3D Systems, Inc. (1988). StereoLithography Interface Specification. - STLファイル形式の公式仕様書、ASCII/Binary STL構造の定義
7. numpy-stl Documentation. (2024). Python library for working with STL files. https://numpy-stl.readthedocs.io/ - STLファイル読込・体積計算のためのPythonライブラリ
8. trimesh Documentation. (2024). Python library for loading and using triangular meshes. https://trimsh.org/ - メッシュ修復・ブーリアン演算・品質評価の包括的ライブラリ

使用ツールとライブラリ

• NumPy (v1.24+): 数値計算ライブラリ - https://numpy.org/
• numpy-stl (v3.0+): STLファイル処理ライブラリ - https://numpy-stl.readthedocs.io/
• trimesh (v4.0+): 3Dメッシュ処理ライブラリ（修復、検証、ブーリアン演算） - https://trimsh.org/
• Matplotlib (v3.7+): データ可視化ライブラリ - https://matplotlib.org/
• SciPy (v1.10+): 科学技術計算ライブラリ（最適化、補間） - https://scipy.org/