第6章: pycalphadによる相図計算実践

TDBファイルから二元系・三元系相図を計算し、平衡相組成・駆動力を求める実践的ワークフローを習得します

📖 読了時間: 20-25 minutes 📊 難易度: 中級〜上級 💻 コード例: 7個 📝 演習問題: 0問

学習目標

この最終章を完了することで、以下の実践的スキルを習得できます：

1. pycalphadの概要と環境構築

1.1 pycalphadとは

pycalphad は、CALPHAD法（CALculation of PHAse Diagrams）に基づいた相図計算を Pythonで実行するためのオープンソースライブラリです。CALPHADデータベース（TDBファイル）を読み込み、 多成分系の相平衡、熱力学量、駆動力などを計算できます。

pycalphadの主要機能

🎯 CALPHAD法の重要性

CALPHAD法は、実験データと第一原理計算を統合して多成分系の熱力学を記述する手法です。 材料開発では、実験で全ての組成・温度を網羅することは不可能ですが、CALPHAD法により 広範な条件での相平衡を予測できます。pycalphadはこの強力な手法をPythonで実装し、 データ駆動型材料設計の基盤となっています。

1.2 インストールとセットアップ

環境構築手順

# ターミナル/コマンドプロンプトで実行

# 基本的なインストール（condaを推奨）
conda install -c conda-forge pycalphad

# または pip でのインストール
pip install pycalphad

# 可視化ライブラリも併せてインストール
pip install matplotlib numpy scipy xarray

# インストール確認
python -c "import pycalphad; print(pycalphad.__version__)"
# 出力例: 0.10.3

# インタラクティブな作業にはJupyter Notebookがおすすめ
pip install jupyter
jupyter notebook

1.3 TDBファイルの取得方法

pycalphadはTDB（ThermoCalc DataBase）形式の熱力学データベースを使用します。 TDBファイルには相の熱力学的記述（ギブスエネルギー関数）が含まれています。

TDBファイルの入手先

⚠️ TDBファイルの制約

商用データベースは再配布禁止のため、本章では公開データベース または デモ用の簡易TDB を使用します。実際の研究では、精度の高い商用データベースの ライセンス取得を推奨します（多くの大学・研究機関で利用可能）。

```mermaid
graph TD
    A[TDBファイル取得] --> B[pycalphadで読み込み]
    B --> C[相図計算]
    B --> D[平衡計算]
    B --> E[駆動力計算]
    C --> F[可視化]
    D --> F
    E --> F
    F --> G[材料設計への応用]

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

2. pycalphadの基本操作

2.1 データベースの読み込みと基本情報の取得

コード例1：TDBファイルの読み込みと基本操作

簡易的なAl-Cu二元系TDBファイルを使った基本操作を示します：

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

# デモ用の簡易Al-Cu TDBファイル（文字列として定義）
# 実際の研究では外部ファイルから読み込む
demo_tdb = """
$ Al-Cu binary system (simplified demo)
ELEMENT AL FCC_A1 26.98 4577.3 28.3 !
ELEMENT CU FCC_A1 63.546 5004.1 33.15 !

FUNCTION GHSERAL 298.15 -7976.15+137.093038*T-24.3671976*T*LN(T)
    -.001884662*T**2-8.77664E-07*T**3+74092*T**(-1); 933.47 Y
    -11276.24+223.048446*T-38.5844296*T*LN(T)+.018531982*T**2
    -5.764227E-06*T**3+74092*T**(-1); 2900 N !

FUNCTION GHSERCU 298.15 -7770.458+130.485235*T-24.112392*T*LN(T)
    -.00265684*T**2+1.29223E-07*T**3+52478*T**(-1); 1358 Y
    -13542.026+183.803828*T-31.38*T*LN(T)+3.64167E+29*T**(-9); 3200 N !

PHASE FCC_A1 % 1 1 !
CONSTITUENT FCC_A1 : AL,CU : !
PARAMETER G(FCC_A1,AL;0) 298.15 +GHSERAL; 6000 N !
PARAMETER G(FCC_A1,CU;0) 298.15 +GHSERCU; 6000 N !
PARAMETER G(FCC_A1,AL,CU;0) 298.15 -53520+7.2*T; 6000 N !

PHASE LIQUID % 1 1 !
CONSTITUENT LIQUID : AL,CU : !
PARAMETER G(LIQUID,AL;0) 298.15 +GHSERAL+11005-11.841*T; 6000 N !
PARAMETER G(LIQUID,CU;0) 298.15 +GHSERCU+12964-9.511*T; 6000 N !
PARAMETER G(LIQUID,AL,CU;0) 298.15 -66200+40*T; 6000 N !
"""

# TDBファイルの読み込み
# 実際のファイルの場合: db = Database('alcu.tdb')
from io import StringIO
db = Database(StringIO(demo_tdb))

print("=== データベース基本情報 ===\n")

# 含まれる元素
print(f"元素: {sorted(db.elements)}")

# 定義されている相
print(f"相: {sorted(db.phases.keys())}")

# 各相の構成情報
print("\n=== 相の詳細情報 ===")
for phase_name in sorted(db.phases.keys()):
    phase = db.phases[phase_name]
    print(f"\n{phase_name}相:")
    print(f"  サブラティス数: {len(phase.constituents)}")
    print(f"  構成元素: {phase.constituents}")

# 元素の基本情報
print("\n=== 元素情報 ===")
for element in sorted(db.elements):
    if element != 'VA':  # 空孔を除く
        print(f"{element}:")
        # 元素のエンタルピー参照状態
        # 実際のTDBにはより詳細な情報が含まれる

print("\n=== pycalphad 変数の確認 ===")
print(f"温度変数: {v.T}")
print(f"圧力変数: {v.P}")
print(f"組成変数（例）: {v.X('AL')}")

print("\nデータベースの読み込みに成功しました。")
print("次は相図計算に進みます。")

2.2 変数と条件の設定

pycalphadでは、variablesモジュールを使って計算条件を設定します。 主要な変数は以下の通りです：

変数記号説明単位
v.TT温度K（ケルビン）
v.PP圧力Pa（パスカル）
v.X('ELEMENT')X_ELEMENT元素のモル分率無次元（0〜1）
v.NN全モル数mol
v.GEGE過剰ギブスエネルギーJ/mol

3. 二元系相図の計算

3.1 Al-Cu二元系相図

コード例2：Al-Cu二元系相図の計算と可視化

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

# 前のコード例で定義したdemo_tdbを使用
from io import StringIO
db = Database(StringIO(demo_tdb))

# 計算条件の設定
components = ['AL', 'CU', 'VA']  # VA（空孔）は必須
phases = list(db.phases.keys())  # ['FCC_A1', 'LIQUID']

# 温度範囲：600 K 〜 1400 K（Al融点付近をカバー）
temperatures = np.linspace(600, 1400, 50)

# 組成範囲：Cu 0 〜 100 at.%（Alの残りを自動計算）
cu_compositions = np.linspace(0, 1, 50)

print("=== Al-Cu二元系相図の計算 ===\n")
print(f"温度範囲: {temperatures[0]:.0f} - {temperatures[-1]:.0f} K")
print(f"組成範囲: Cu {cu_compositions[0]:.1%} - {cu_compositions[-1]:.1%}")
print(f"計算点数: {len(temperatures)} × {len(cu_compositions)} = {len(temperatures)*len(cu_compositions)}")

# 平衡計算の実行
try:
    eq_result = equilibrium(
        db,
        components,
        phases,
        {
            v.X('CU'): cu_compositions,
            v.T: temperatures,
            v.P: 101325,  # 標準圧力（Pa）
            v.N: 1.0      # 全モル数（任意の値でOK）
        },
        output='GM'  # モルギブスエネルギーを出力
    )

    print("\n計算完了！")

    # 結果の可視化
    fig, ax = plt.subplots(figsize=(10, 7))

    # 相境界の抽出と描画
    # eq_result.Phase には各点での安定相情報が格納されている
    # NP（安定相の数）が変化する点が相境界

    # 簡易的な可視化（実際にはより洗練された方法を使用）
    for phase_name in phases:
        # 各相が安定な領域を色分け表示
        phase_fraction = eq_result.NP.where(
            eq_result.Phase == phase_name
        ).values

        # 相図の描画（カラーマップ）
        im = ax.contourf(
            cu_compositions * 100,  # at.% に変換
            temperatures,
            phase_fraction.T,
            levels=[0, 0.01, 1],
            colors=['white', 'lightblue'] if phase_name == 'LIQUID' else ['white', 'lightcoral'],
            alpha=0.5
        )

    ax.set_xlabel('Cu (at.%)', fontsize=12, fontweight='bold')
    ax.set_ylabel('Temperature (K)', fontsize=12, fontweight='bold')
    ax.set_title('Al-Cu Binary Phase Diagram (Simplified Demo)',
                 fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)

    # 凡例を追加
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor='lightcoral', alpha=0.5, label='FCC_A1 (Solid)'),
        Patch(facecolor='lightblue', alpha=0.5, label='LIQUID')
    ]
    ax.legend(handles=legend_elements, loc='upper right', fontsize=10)

    plt.tight_layout()
    plt.savefig('alcu_phase_diagram_demo.png', dpi=150, bbox_inches='tight')
    plt.show()

    print("\n相図を保存しました: alcu_phase_diagram_demo.png")

except Exception as e:
    print(f"\n計算中にエラーが発生しました: {e}")
    print("デモ用簡易TDBでは完全な相図計算ができない場合があります。")
    print("実際の研究では高精度なTDBファイルを使用してください。")

# 実際の研究でのベストプラクティス
print("\n=== 実際の研究での推奨事項 ===")
print("1. 高精度な商用TDB（TCAL等）またはpycalphad-dataを使用")
print("2. 計算点数を増やして相境界を滑らかに（100×100以上）")
print("3. binplot()関数を使った専用の相図描画")
print("4. 実験データと比較して精度を検証")

3.2 Fe-C二元系相図（鉄鋼材料の基礎）

コード例3：Fe-C二元系相図の計算

鉄鋼材料で最も重要なFe-C系の相図を計算します：

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

# Fe-C系の簡易TDB（実際にはTCAL5等の精密データベースを使用）
# ここではデモ用の概念的なコードを示します

print("=== Fe-C二元系相図の計算例 ===\n")

# 実際の計算には適切なTDBファイルが必要
# 例: db = Database('tcfe9.tdb')  # Thermo-Calc鉄鋼データベース

# デモモード：計算手順の説明
print("【計算手順】")
print("1. Fe-C系TDBファイルの読み込み")
print("   db = Database('tcfe9.tdb')")
print()
print("2. 計算条件の設定")
print("   components = ['FE', 'C', 'VA']")
print("   phases = ['LIQUID', 'BCC_A2', 'FCC_A1', 'CEMENTITE']")
print("   # BCC_A2: フェライト（α鉄）")
print("   # FCC_A1: オーステナイト（γ鉄）")
print("   # CEMENTITE: セメンタイト（Fe3C）")
print()
print("3. 温度・組成範囲の設定")
print("   temperatures = np.linspace(800, 1800, 100)  # K")
print("   carbon_content = np.linspace(0, 0.08, 100)  # C濃度 0-8 wt.%")
print()
print("4. 平衡計算")
print("   eq = equilibrium(db, components, phases, {")
print("       v.X('C'): carbon_content,")
print("       v.T: temperatures,")
print("       v.P: 101325")
print("   })")
print()
print("5. 相図の可視化")
print("   # 液相線、固相線、共析温度（727°C）などが可視化される")
print()

# 重要な相変態温度の理論値
print("【Fe-C系の重要な相変態温度】")
phase_transformations = {
    "共析温度（Eutectoid）": "727°C (1000 K)",
    "共析組成": "0.76 wt.% C",
    "包析温度（Peritectic）": "1493°C (1766 K)",
    "共晶温度（Eutectic）": "1147°C (1420 K)",
    "α-γ変態温度（純Fe）": "910°C (1183 K)"
}

for name, value in phase_transformations.items():
    print(f"  {name}: {value}")

# 実際の相図の特徴を可視化（模式図）
fig, ax = plt.subplots(figsize=(11, 8))

# 模式的な相図の描画（実データではなく概念図）
# 実際の計算結果を表示するには上記のequilibrium()が必要

# 温度軸（縦軸）
T_liquid = np.array([1768, 1493, 1147])  # K
T_alpha = np.array([1183, 1000, 1000])
T_gamma = np.array([1493, 1147, 1000])

# 組成軸（横軸、wt.% C）
C_liquid = np.array([0, 0.17, 4.3])
C_alpha = np.array([0, 0.02, 0.02])
C_gamma = np.array([0.17, 2.14, 0.76])

ax.plot(C_liquid, T_liquid, 'b-', linewidth=2, label='Liquidus')
ax.plot(C_alpha, T_alpha, 'r-', linewidth=2, label='α (BCC) boundary')
ax.plot(C_gamma, T_gamma, 'g-', linewidth=2, label='γ (FCC) boundary')

# 重要な点をマーク
ax.plot(0.76, 1000, 'ko', markersize=10, label='Eutectoid point')
ax.plot(0.17, 1766, 'ms', markersize=10, label='Peritectic point')
ax.plot(4.3, 1420, 'c^', markersize=10, label='Eutectic point')

# 相領域のラベル
ax.text(0.1, 1400, 'LIQUID', fontsize=12, fontweight='bold', color='blue')
ax.text(0.4, 1300, 'L + γ', fontsize=11, color='green')
ax.text(0.01, 1050, 'α', fontsize=12, fontweight='bold', color='red')
ax.text(0.3, 1100, 'γ', fontsize=12, fontweight='bold', color='green')
ax.text(0.4, 950, 'α + Fe₃C', fontsize=11, color='purple')
ax.text(1.5, 950, 'γ + Fe₃C', fontsize=11, color='purple')

ax.set_xlabel('Carbon content (wt.%)', fontsize=13, fontweight='bold')
ax.set_ylabel('Temperature (K)', fontsize=13, fontweight='bold')
ax.set_title('Fe-C Binary Phase Diagram (Schematic)',
             fontsize=15, fontweight='bold')
ax.set_xlim(0, 2.5)
ax.set_ylim(900, 1900)
ax.grid(True, alpha=0.3)
ax.legend(loc='upper right', fontsize=10)

plt.tight_layout()
plt.savefig('fec_phase_diagram_schematic.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n模式図を保存しました: fec_phase_diagram_schematic.png")
print("\n【注意】上記は模式図です。実際の研究ではTDBファイルからの計算が必要です。")

Fe-C相図の材料科学的意義

Fe-C二元系相図は、鉄鋼材料の熱処理（焼入れ、焼戻し、焼なまし等）の基礎となります。 例えば、共析温度（727°C、0.76 wt.% C） を境に、 オーステナイト（γ相）がフェライト（α相）とセメンタイト（Fe₃C）の混合組織である パーライト に変態します。この知識は、鋼の機械的性質を制御する上で不可欠です。

4. 三元系相図の計算

4.1 等温断面図（Isothermal Section）

三元系では、温度を固定した等温断面図 を計算することが一般的です。 三角図（Gibbs三角形）上に相領域を表示します。

コード例4：Al-Cu-Mg三元系の等温断面図

アルミニウム合金で重要なAl-Cu-Mg系の相図を計算します：

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

print("=== Al-Cu-Mg三元系等温断面図の計算例 ===\n")

# 実際の計算には三元系TDBファイルが必要
# 例: db = Database('alzn_mey.tdb')  # Al-Zn-Mg系などのTDB

# 計算手順の説明（デモモード）
print("【計算手順】")
print("1. 三元系TDBファイルの読み込み")
print("   db = Database('al-cu-mg.tdb')")
print()
print("2. 等温断面の温度を設定（例：600 K）")
print("   temperature = 600  # K (約327°C)")
print()
print("3. 三角図上の組成グリッドを生成")
print("   # Al-Cu-Mg三角図上の点を均等に配置")
print("   # X(Al) + X(Cu) + X(Mg) = 1 の制約")
print()
print("4. 各組成点で平衡計算")
print("   eq = equilibrium(db, components, phases, {")
print("       v.X('CU'): cu_fractions,")
print("       v.X('MG'): mg_fractions,")
print("       v.T: temperature,")
print("       v.P: 101325")
print("   })")
print()
print("5. 三角図上に相領域を可視化")
print()

# 三角図の座標変換関数（Gibbs三角形用）
def ternary_to_cartesian(a, b, c):
    """
    三成分（a, b, c）の組成を三角図のデカルト座標（x, y）に変換

    Parameters:
    -----------
    a, b, c : float or array
        3成分の組成（a + b + c = 1）

    Returns:
    --------
    x, y : float or array
        三角図上の座標
    """
    x = 0.5 * (2 * b + c) / (a + b + c)
    y = (np.sqrt(3) / 2) * c / (a + b + c)
    return x, y

# 三角図の枠を描画
fig, ax = plt.subplots(figsize=(10, 9))

# 三角形の頂点（Al, Cu, Mg）
vertices = np.array([
    [0, 0],           # Al (100%)
    [1, 0],           # Cu (100%)
    [0.5, np.sqrt(3)/2]  # Mg (100%)
])

triangle = plt.Polygon(vertices, fill=False, edgecolor='black', linewidth=2)
ax.add_patch(triangle)

# 頂点のラベル
ax.text(0, -0.05, 'Al', fontsize=14, fontweight='bold', ha='center', va='top')
ax.text(1, -0.05, 'Cu', fontsize=14, fontweight='bold', ha='center', va='top')
ax.text(0.5, np.sqrt(3)/2 + 0.05, 'Mg', fontsize=14, fontweight='bold',
        ha='center', va='bottom')

# グリッド線（10%刻み）
for i in range(1, 10):
    frac = i / 10

    # Al-Cu軸に平行な線（Mg一定）
    x1, y1 = ternary_to_cartesian(1-frac, 0, frac)
    x2, y2 = ternary_to_cartesian(0, 1-frac, frac)
    ax.plot([x1, x2], [y1, y2], 'k-', linewidth=0.5, alpha=0.3)

    # Cu-Mg軸に平行な線（Al一定）
    x1, y1 = ternary_to_cartesian(frac, 1-frac, 0)
    x2, y2 = ternary_to_cartesian(frac, 0, 1-frac)
    ax.plot([x1, x2], [y1, y2], 'k-', linewidth=0.5, alpha=0.3)

    # Mg-Al軸に平行な線（Cu一定）
    x1, y1 = ternary_to_cartesian(1-frac, frac, 0)
    x2, y2 = ternary_to_cartesian(0, frac, 1-frac)
    ax.plot([x1, x2], [y1, y2], 'k-', linewidth=0.5, alpha=0.3)

# 模式的な相領域を表示（実際の計算結果の代わり）
# 実際にはequilibrium()の結果からcontour plotを作成

# FCC_A1相領域（Al-richな固溶体）
al_rich_x = [0, 0.3, 0.15, 0]
al_rich_y = [0, 0, 0.13, 0]
ax.fill(al_rich_x, al_rich_y, color='lightcoral', alpha=0.5,
        label='FCC_A1 (Al-rich)')

# その他の相領域（概念的）
ax.text(0.15, 0.05, 'α-Al', fontsize=11, fontweight='bold', ha='center')
ax.text(0.6, 0.3, 'Intermetallic\nphases', fontsize=10, ha='center', style='italic')

ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.1, 1.0)
ax.set_aspect('equal')
ax.axis('off')
ax.set_title('Al-Cu-Mg Ternary Isothermal Section at 600 K (Schematic)',
             fontsize=14, fontweight='bold', pad=20)

# 凡例
ax.legend(loc='upper right', fontsize=10)

plt.tight_layout()
plt.savefig('al-cu-mg_ternary_schematic.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n三角図を保存しました: al-cu-mg_ternary_schematic.png")
print("\n【注意】上記は模式図です。実際の研究ではTDBファイルからの計算が必要です。")

print("\n=== 三元系計算の実践的ポイント ===")
print("1. 計算時間：三元系は二元系より計算コストが高い（数分〜数時間）")
print("2. グリッド解像度：50×50程度が標準、高精度なら100×100")
print("3. 相の選択：計算対象の相を適切に選ぶ（全相だと重くなる）")
print("4. 可視化：三角図プロット、等高線図、相分率マップ等")

4.2 垂直断面図（Vertical Section）

三元系の特定の組成比を固定し、温度を変化させた垂直断面図も重要です。 例えば、Al-Cu-Mg系でMg濃度を固定し、Al-Cu比と温度の関係を調べることができます。

5. 平衡相組成と相分率の計算

5.1 指定条件での平衡計算

コード例5：平衡相組成と相分率の計算

特定の温度・組成での平衡状態を詳細に解析します：

import pycalphad as pyc
from pycalphad import Database, equilibrium, variables as v
import numpy as np
import pandas as pd

print("=== 平衡相組成と相分率の計算 ===\n")

# 前のコード例のdemo_tdb（Al-Cu系）を使用
from io import StringIO
db = Database(StringIO(demo_tdb))

components = ['AL', 'CU', 'VA']
phases = list(db.phases.keys())

# 計算条件：特定の温度と組成
temperature = 900  # K （約627°C）
cu_content = 0.3   # Cu 30 at.%

print(f"計算条件:")
print(f"  温度: {temperature} K ({temperature - 273.15:.1f}°C)")
print(f"  組成: Al-{cu_content*100:.0f}at.%Cu")
print(f"  圧力: 101325 Pa (1 atm)")
print()

try:
    # 平衡計算
    eq = equilibrium(
        db, components, phases,
        {
            v.X('CU'): cu_content,
            v.T: temperature,
            v.P: 101325,
            v.N: 1.0
        },
        output=['GM', 'NP', 'X']  # モルギブスエネルギー、相分率、組成
    )

    print("平衡計算完了！\n")
    print("=== 計算結果 ===\n")

    # 安定相の抽出
    stable_phases = []
    for phase_name in phases:
        # 相分率（NP）が0より大きい相が安定相
        phase_fraction = float(eq.NP.sel(Phase=phase_name).values)
        if phase_fraction > 1e-6:  # 微小な数値誤差を除外
            stable_phases.append({
                'Phase': phase_name,
                'Fraction': phase_fraction
            })

    if stable_phases:
        print("安定相:")
        for phase_data in stable_phases:
            print(f"  {phase_data['Phase']}: {phase_data['Fraction']:.4f} ({phase_data['Fraction']*100:.2f}%)")
    else:
        print("安定相が見つかりませんでした（計算エラーの可能性）")

    print()

    # 各相の組成
    print("各相の組成:")
    for phase_data in stable_phases:
        phase_name = phase_data['Phase']
        print(f"\n  {phase_name}相:")

        # Al, Cuそれぞれの組成
        try:
            al_content = float(eq.X.sel(Phase=phase_name, component='AL').values)
            cu_content_phase = float(eq.X.sel(Phase=phase_name, component='CU').values)

            print(f"    Al: {al_content:.4f} ({al_content*100:.2f} at.%)")
            print(f"    Cu: {cu_content_phase:.4f} ({cu_content_phase*100:.2f} at.%)")
        except:
            print("    組成データの取得に失敗")

    # ギブスエネルギー
    print(f"\nシステム全体のモルギブスエネルギー:")
    gm_total = float(eq.GM.values)
    print(f"  {gm_total:.2f} J/mol")

    print("\n=== 実際の材料設計への応用 ===")
    print("この情報から以下が予測できます：")
    print("1. どの相が析出するか（相分率から判断）")
    print("2. 各相の組成（偏析の程度）")
    print("3. 熱処理後の組織（相の安定性）")
    print("4. 機械的性質への影響（相の種類と分率）")

except Exception as e:
    print(f"計算エラー: {e}")
    print("デモ用TDBでは詳細な計算ができない場合があります。")

# 温度を変えて相分率の変化を追跡
print("\n\n=== 温度を変化させた相分率の変化 ===\n")

temperatures_range = np.linspace(700, 1100, 20)
results = []

for T in temperatures_range:
    try:
        eq_temp = equilibrium(
            db, components, phases,
            {v.X('CU'): 0.3, v.T: T, v.P: 101325, v.N: 1.0},
            output='NP'
        )

        fcc_fraction = float(eq_temp.NP.sel(Phase='FCC_A1').values)
        liquid_fraction = float(eq_temp.NP.sel(Phase='LIQUID').values)

        results.append({
            'T (K)': T,
            'T (°C)': T - 273.15,
            'FCC_A1': fcc_fraction,
            'LIQUID': liquid_fraction
        })
    except:
        pass

if results:
    df = pd.DataFrame(results)
    print(df.to_string(index=False))

    # 可視化
    import matplotlib.pyplot as plt

    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(df['T (K)'], df['FCC_A1'], 'ro-', linewidth=2,
            markersize=6, label='FCC_A1 (Solid)')
    ax.plot(df['T (K)'], df['LIQUID'], 'bo-', linewidth=2,
            markersize=6, label='LIQUID')

    ax.set_xlabel('Temperature (K)', fontsize=12, fontweight='bold')
    ax.set_ylabel('Phase Fraction', fontsize=12, fontweight='bold')
    ax.set_title('Phase Fraction vs Temperature (Al-30at.%Cu)',
                 fontsize=14, fontweight='bold')
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)
    ax.set_ylim(-0.05, 1.05)

    plt.tight_layout()
    plt.savefig('phase_fraction_vs_temperature.png', dpi=150, bbox_inches='tight')
    plt.show()

    print("\n相分率グラフを保存しました: phase_fraction_vs_temperature.png")

5.2 レバー則による相分率計算

二相共存領域では、レバー則（lever rule） を使って各相の量を求めることができます。 pycalphadは自動的にこれを計算しますが、原理を理解することが重要です。

🎯 レバー則の原理

二相α, βが共存する場合、全体組成$X_0$における各相の分率は：

$f_\alpha = \frac{X_\beta - X_0}{X_\beta - X_\alpha}, \quad f_\beta = \frac{X_0 - X_\alpha}{X_\beta - X_\alpha}$

ここで、$X_\alpha$, $X_\beta$はそれぞれα相、β相の組成、$X_0$は全体組成です。 これは「てこの原理」と同じ数学的関係です。

6. 駆動力（Driving Force）の計算

6.1 駆動力とは

駆動力（driving force） は、ある相から別の相への変態がどれだけ進行しやすいかを 示す熱力学的指標です。化学ポテンシャルの差として定義されます。

駆動力の定義

相αから相βへの変態駆動力は：

$\Delta G = G_\beta - G_\alpha$

駆動力が大きいほど、相変態の速度が速くなります（ただし、速度論的障壁も考慮が必要）。

コード例6：相変態駆動力の計算

Al-Cu系でFCC相からLIQUID相への融解駆動力を計算します：

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

print("=== 相変態駆動力の計算 ===\n")

# demo_tdbを使用
from io import StringIO
db = Database(StringIO(demo_tdb))

components = ['AL', 'CU', 'VA']
phases = list(db.phases.keys())

# 計算条件
cu_composition = 0.2  # Cu 20 at.%
temperatures = np.linspace(700, 1200, 30)

print(f"組成: Al-{cu_composition*100:.0f}at.%Cu")
print(f"温度範囲: {temperatures[0]:.0f} - {temperatures[-1]:.0f} K")
print()

# 各温度での各相のギブスエネルギーを計算
fcc_energies = []
liquid_energies = []
driving_forces = []

for T in temperatures:
    try:
        # FCC相のギブスエネルギー
        eq_fcc = equilibrium(
            db, components, ['FCC_A1'],  # FCC相のみ
            {v.X('CU'): cu_composition, v.T: T, v.P: 101325},
            output='GM'
        )
        G_fcc = float(eq_fcc.GM.values)
        fcc_energies.append(G_fcc)

        # LIQUID相のギブスエネルギー
        eq_liquid = equilibrium(
            db, components, ['LIQUID'],  # LIQUID相のみ
            {v.X('CU'): cu_composition, v.T: T, v.P: 101325},
            output='GM'
        )
        G_liquid = float(eq_liquid.GM.values)
        liquid_energies.append(G_liquid)

        # 駆動力：LIQUID - FCC（正なら液相が安定）
        driving_force = G_liquid - G_fcc
        driving_forces.append(driving_force)

    except Exception as e:
        fcc_energies.append(np.nan)
        liquid_energies.append(np.nan)
        driving_forces.append(np.nan)

print("計算完了！\n")

# 結果の可視化
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

# サブプロット1：ギブスエネルギー
ax1.plot(temperatures, fcc_energies, 'r-', linewidth=2.5, label='FCC_A1 (Solid)')
ax1.plot(temperatures, liquid_energies, 'b-', linewidth=2.5, label='LIQUID')
ax1.axhline(y=0, color='gray', linestyle='--', linewidth=1, alpha=0.5)
ax1.set_xlabel('Temperature (K)', fontsize=12, fontweight='bold')
ax1.set_ylabel('Gibbs Energy (J/mol)', fontsize=12, fontweight='bold')
ax1.set_title(f'Gibbs Energy vs Temperature (Al-{cu_composition*100:.0f}at.%Cu)',
              fontsize=13, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# サブプロット2：駆動力
ax2.plot(temperatures, driving_forces, 'g-', linewidth=2.5, label='ΔG (LIQUID - FCC)')
ax2.axhline(y=0, color='red', linestyle='--', linewidth=2, alpha=0.7, label='Equilibrium')
ax2.fill_between(temperatures, 0, driving_forces,
                 where=np.array(driving_forces) < 0, alpha=0.3, color='blue',
                 label='FCC stable region')
ax2.fill_between(temperatures, 0, driving_forces,
                 where=np.array(driving_forces) > 0, alpha=0.3, color='red',
                 label='LIQUID stable region')
ax2.set_xlabel('Temperature (K)', fontsize=12, fontweight='bold')
ax2.set_ylabel('Driving Force (J/mol)', fontsize=12, fontweight='bold')
ax2.set_title(f'Driving Force for Melting (Al-{cu_composition*100:.0f}at.%Cu)',
              fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('driving_force_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

print("駆動力解析グラフを保存しました: driving_force_analysis.png\n")

# 融点の推定（駆動力がゼロになる温度）
try:
    # ゼロクロス点を探す
    for i in range(len(driving_forces) - 1):
        if not np.isnan(driving_forces[i]) and not np.isnan(driving_forces[i+1]):
            if driving_forces[i] * driving_forces[i+1] < 0:  # 符号が変わる
                T_melting = (temperatures[i] + temperatures[i+1]) / 2
                print(f"推定融点: {T_melting:.1f} K ({T_melting - 273.15:.1f}°C)\n")
                break
except:
    pass

print("=== 駆動力計算の応用 ===")
print("1. 相変態の熱力学的実行可能性の評価")
print("2. 過冷却・過加熱の程度の定量化")
print("3. 核生成・成長速度の推定（速度論と組み合わせ）")
print("4. 熱処理条件の最適化")
print("5. 準安定相の安定性評価")

6.2 駆動力と相変態速度の関係

駆動力は相変態の「熱力学的な推進力」を表しますが、実際の変態速度は 速度論 （拡散、界面移動など）にも依存します。 両者を組み合わせた解析が、実際の材料プロセス設計には不可欠です。

7. 実践的ワークフローと応用例

7.1 完全な研究ワークフロー

```mermaid
graph TD
    A[研究目的の明確化] --> B[適切なTDB選択]
    B --> C[pycalphadで計算]
    C --> D{計算対象}
    D --> E[相図計算]
    D --> F[平衡計算]
    D --> G[駆動力計算]
    E --> H[結果の可視化]
    F --> H
    G --> H
    H --> I[実験データと比較]
    I --> J{精度OK?}
    J -->|Yes| K[材料設計への応用]
    J -->|No| L[TDB改善/パラメータ調整]
    L --> C
    K --> M[実験での検証]
    M --> N[論文・特許]

    style A fill:#e3f2fd
    style C fill:#fff3e0
    style H fill:#e8f5e9
    style K fill:#fce4ec
    style N fill:#fff9c4
```

コード例7：実践的ワークフロー統合スクリプト

TDB取得から可視化、材料設計までの完全なワークフローを示します：

import pycalphad as pyc
from pycalphad import Database, equilibrium, variables as v
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from datetime import datetime

print("=" * 60)
print("pycalphad 実践的ワークフロー統合スクリプト")
print("=" * 60)
print()

# ============================================
# ステップ1: プロジェクト設定
# ============================================
print("【ステップ1】プロジェクト設定\n")

project_name = "Al-Cu Alloy Design"
target_composition = "Al-20at.%Cu"
target_temperature_range = (600, 1200)  # K
objective = "最適熱処理条件の決定"

print(f"プロジェクト名: {project_name}")
print(f"対象合金: {target_composition}")
print(f"温度範囲: {target_temperature_range[0]} - {target_temperature_range[1]} K")
print(f"目的: {objective}")
print()

# ============================================
# ステップ2: データベース準備
# ============================================
print("【ステップ2】TDBファイルの読み込み\n")

# 実際の研究では外部TDBファイルを使用
# db = Database('path/to/your/database.tdb')

# デモ用
from io import StringIO
db = Database(StringIO(demo_tdb))

components = ['AL', 'CU', 'VA']
phases = list(db.phases.keys())

print(f"使用するTDB: demo_al_cu.tdb")
print(f"元素: {[c for c in components if c != 'VA']}")
print(f"相: {phases}")
print()

# ============================================
# ステップ3: 相図計算
# ============================================
print("【ステップ3】相図計算\n")

T_range = np.linspace(target_temperature_range[0],
                      target_temperature_range[1], 40)
X_Cu_range = np.linspace(0, 0.5, 40)  # 0-50 at.% Cu

print(f"計算中... (温度{len(T_range)}点 × 組成{len(X_Cu_range)}点)")

try:
    eq_phase_diagram = equilibrium(
        db, components, phases,
        {v.X('CU'): X_Cu_range, v.T: T_range, v.P: 101325},
        output='GM'
    )
    print("✓ 相図計算完了")
except Exception as e:
    print(f"✗ 相図計算エラー: {e}")
    eq_phase_diagram = None

print()

# ============================================
# ステップ4: 特定条件での平衡計算
# ============================================
print("【ステップ4】特定条件での平衡計算\n")

specific_T = 900  # K
specific_X_Cu = 0.2  # 20 at.%

print(f"条件: {specific_T} K, Cu {specific_X_Cu*100:.0f} at.%")

try:
    eq_specific = equilibrium(
        db, components, phases,
        {v.X('CU'): specific_X_Cu, v.T: specific_T, v.P: 101325},
        output=['GM', 'NP', 'X']
    )

    print("\n安定相:")
    for phase_name in phases:
        fraction = float(eq_specific.NP.sel(Phase=phase_name).values)
        if fraction > 1e-6:
            print(f"  {phase_name}: {fraction:.4f} ({fraction*100:.2f}%)")

    print("✓ 平衡計算完了")
except Exception as e:
    print(f"✗ 平衡計算エラー: {e}")
    eq_specific = None

print()

# ============================================
# ステップ5: 駆動力解析
# ============================================
print("【ステップ5】駆動力解析\n")

T_df_range = np.linspace(700, 1100, 25)
driving_forces_analysis = []

for T in T_df_range:
    try:
        eq_fcc = equilibrium(db, components, ['FCC_A1'],
                            {v.X('CU'): 0.2, v.T: T, v.P: 101325}, output='GM')
        eq_liq = equilibrium(db, components, ['LIQUID'],
                            {v.X('CU'): 0.2, v.T: T, v.P: 101325}, output='GM')

        G_fcc = float(eq_fcc.GM.values)
        G_liq = float(eq_liq.GM.values)
        df_value = G_liq - G_fcc

        driving_forces_analysis.append({
            'T (K)': T,
            'ΔG (J/mol)': df_value
        })
    except:
        pass

if driving_forces_analysis:
    df_analysis = pd.DataFrame(driving_forces_analysis)
    print(df_analysis.head(10).to_string(index=False))
    print("...")
    print("✓ 駆動力解析完了")
else:
    print("✗ 駆動力解析エラー")

print()

# ============================================
# ステップ6: 結果の可視化
# ============================================
print("【ステップ6】結果の可視化\n")

fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# (1) 相図（簡易版）
ax1 = axes[0, 0]
if eq_phase_diagram is not None:
    # 相図の描画（実際にはより詳細な処理が必要）
    ax1.text(0.5, 0.5, 'Phase Diagram\n(Requires detailed plotting)',
             ha='center', va='center', fontsize=12, transform=ax1.transAxes)
ax1.set_title('Phase Diagram', fontweight='bold')
ax1.set_xlabel('Cu (at.%)')
ax1.set_ylabel('Temperature (K)')
ax1.grid(True, alpha=0.3)

# (2) 相分率 vs 温度
ax2 = axes[0, 1]
if not driving_forces_analysis:
    ax2.text(0.5, 0.5, 'No data', ha='center', va='center',
             transform=ax2.transAxes)
else:
    # ダミーデータ
    ax2.plot([700, 900, 1100], [1, 0.7, 0], 'ro-', label='FCC', linewidth=2)
    ax2.plot([700, 900, 1100], [0, 0.3, 1], 'bo-', label='LIQUID', linewidth=2)
ax2.set_title('Phase Fraction vs Temperature', fontweight='bold')
ax2.set_xlabel('Temperature (K)')
ax2.set_ylabel('Phase Fraction')
ax2.legend()
ax2.grid(True, alpha=0.3)

# (3) 駆動力
ax3 = axes[1, 0]
if driving_forces_analysis:
    ax3.plot(df_analysis['T (K)'], df_analysis['ΔG (J/mol)'],
             'g-', linewidth=2.5)
    ax3.axhline(y=0, color='red', linestyle='--', linewidth=2)
    ax3.fill_between(df_analysis['T (K)'], 0, df_analysis['ΔG (J/mol)'],
                     alpha=0.3)
ax3.set_title('Driving Force (LIQUID - FCC)', fontweight='bold')
ax3.set_xlabel('Temperature (K)')
ax3.set_ylabel('ΔG (J/mol)')
ax3.grid(True, alpha=0.3)

# (4) 材料設計指針
ax4 = axes[1, 1]
ax4.axis('off')
design_guidelines = f"""
材料設計指針

【最適熱処理条件】
• 溶体化処理: 950-1000 K
• 時効処理: 450-500 K
• 冷却速度: 急冷推奨

【予想される析出相】
• FCC_A1 (α-Al matrix)
• θ相 (Al₂Cu, 析出強化)

【機械的性質予測】
• 高強度化: θ相析出により
• 延性: マトリックスの性質に依存

【次のステップ】
1. 実験での検証
2. 微細組織観察（SEM/TEM）
3. 機械試験（引張、硬さ）
"""
ax4.text(0.05, 0.95, design_guidelines,
         fontsize=10, verticalalignment='top',
         fontfamily='monospace',
         bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
ax4.set_title('Material Design Guidelines', fontweight='bold', loc='left')

plt.tight_layout()
plt.savefig('workflow_comprehensive_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

print("✓ 総合解析図を保存: workflow_comprehensive_analysis.png")
print()

# ============================================
# ステップ7: レポート生成
# ============================================
print("【ステップ7】レポート生成\n")

report = f"""
{'='*60}
pycalphad 計算レポート
{'='*60}

プロジェクト名: {project_name}
実行日時: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}

【計算条件】
- 対象系: {target_composition}
- 温度範囲: {target_temperature_range[0]} - {target_temperature_range[1]} K
- 圧力: 101325 Pa (1 atm)

【使用データベース】
- TDB: demo_al_cu.tdb
- 相: {', '.join(phases)}

【計算結果サマリー】
- 相図計算: {'成功' if eq_phase_diagram is not None else '失敗'}
- 平衡計算: {'成功' if eq_specific is not None else '失敗'}
- 駆動力解析: {'成功' if driving_forces_analysis else '失敗'}

【材料設計への提言】
{objective}に向けて、以下の条件を推奨します：
1. 溶体化処理: 950-1000 K で完全溶解
2. 急冷: 室温付近の準安定状態を保持
3. 時効処理: 450-500 K で析出強化相を生成

【参考文献】
[1] Pycalphad Documentation: https://pycalphad.org/
[2] CALPHAD methodology review papers
[3] Al-Cu binary phase diagram experimental data

{'='*60}
"""

print(report)

# レポートをファイルに保存
with open('pycalphad_calculation_report.txt', 'w', encoding='utf-8') as f:
    f.write(report)

print("✓ レポートを保存: pycalphad_calculation_report.txt")
print()

print("=" * 60)
print("ワークフロー完了！")
print("=" * 60)

7.2 ベストプラクティス

pycalphad使用時のベストプラクティス

8. 演習問題

演習問題

演習1：基本的な相図計算

公開されているTDBファイル（pycalphad-dataまたはTCAL）を使って、 興味のある二元系（例：Ni-Al, Ti-Al, Cu-Zn等）の相図を計算し、 以下の情報を抽出しなさい：

また、計算した相図を適切に可視化し、論文品質の図を作成しなさい。

演習2：平衡相組成の温度依存性

Al-4wt.%Cu合金を例に、温度を300 K から 900 K まで変化させたときの：

  1. 各相の相分率の変化をプロット
  2. 各相のCu濃度の変化をプロット
  3. 固溶限（solvus）温度を決定
  4. θ相（Al₂Cu）が析出する温度範囲を特定

これらの情報から、最適な時効処理温度を提案しなさい。

演習3：駆動力と過冷却度の関係

純Alの融解を例に、以下を計算しなさい：

  1. 平衡融点（ΔG = 0となる温度）を求める
  2. 融点より10 K、20 K、50 K低い温度での駆動力を計算
  3. 駆動力と過冷却度の関係をプロット
  4. 古典的核生成理論と組み合わせて、臨界核半径を推定

（ヒント：臨界核半径 $r^* = -2\gamma/\Delta G_v$、 界面エネルギー$\gamma \approx 0.1$ J/m²を使用）

演習4：実践プロジェクト

自分の研究テーマに関連する材料系を選び、以下の完全なワークフローを実行しなさい：

  1. 適切なTDBファイルの選定と取得
  2. 相図の計算と文献値との比較
  3. 目標組成での平衡計算（複数温度）
  4. 駆動力解析（必要に応じて）
  5. 結果の総合的可視化（複数のグラフ）
  6. 材料設計への提言をまとめたレポート作成

最終的に、学会発表または論文投稿に使用できるレベルの 図表とレポートを作成することを目標とします。

🎓 シリーズ完了おめでとうございます！

材料熱力学入門シリーズの全6章を完了しました。
あなたは今、pycalphadを使った実践的な相図計算と
熱力学解析のスキルを習得しています。

まとめと次のステップ

本シリーズで習得したスキル

全6章を通じて、以下の知識とスキルを体系的に習得しました：

理論的基盤（第1-3章）

実践的技術（第4-6章）

研究への応用力

次に学ぶべき内容

推奨される発展学習

1. 材料熱力学の高度なトピック

2. 速度論との統合

3. 第一原理計算との連携

4. Materials Informatics（MI）

5. 実験技術との統合

学習の継続のために

推奨リソース

公式ドキュメント

教科書

オンラインコミュニティ

学術論文

最終メッセージ

材料熱力学は、材料科学の最も基盤的な学問領域の一つです。 「なぜこの相が安定なのか」「どの条件でどの相が共存するのか」という 本質的な問いに答えるための強力なツールです。

本シリーズでは、基礎理論からpycalphadを使った実践的計算まで、 体系的に学習しました。特に最終章では、TDBファイルを使った リアルな研究ワークフローを体験しました。

ここで習得した知識とスキルは、合金設計、プロセス最適化、 新材料探索、Materials Informaticsなど、あらゆる材料研究の基盤となります。

これからは、興味のある材料系や研究テーマに対して、 このシリーズで学んだ技術を積極的に応用してください。 pycalphadは強力なツールであり、あなたの研究を大きく加速させる可能性を秘めています。

継続的な学習と実践が、真のスキル習得への道です。
本シリーズが、あなたの材料科学研究の新たなスタートとなることを願っています。

🎓 材料熱力学入門シリーズ 完 🎓

免責事項

← 前の章 シリーズ目次に戻る