第1章: 熱力学の基本法則と熱力学ポテンシャル

import numpy as np import matplotlib.pyplot as plt # 物理定数 R = 8.314 # J/(mol·K) n = 1.0 # mol gamma = 1.4 # 比熱比（空気） # 初期状態 P1 = 1e5 # Pa (1 atm) V1 = 0.0224 # m³ (標準状態) T1 = 273.15 # K # 体積範囲 V_range = np.linspace(0.01, 0.05, 200) # 等温過程 (PV = nRT = const) def isothermal_process(V, n, R, T): """等温過程: PV = nRT""" return n * R * T / V # 断熱過程 (PV^γ = const) def adiabatic_process(V, P1, V1, gamma): """断熱過程: PV^γ = const""" return P1 * (V1 / V)**gamma # 等圧過程 def isobaric_process(V, P): """等圧過程: P = const""" return P * np.ones_like(V) # 計算 P_isothermal = isothermal_process(V_range, n, R, T1) P_adiabatic = adiabatic_process(V_range, P1, V1, gamma) P_isobaric = isobaric_process(V_range, P1) # 可視化 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # PV図 ax1 = axes[0] ax1.plot(V_range * 1000, P_isothermal / 1e5, 'b-', linewidth=2, label=f'等温 (T = {T1:.1f} K)') ax1.plot(V_range * 1000, P_adiabatic / 1e5, 'r-', linewidth=2, label=f'断熱 (γ = {gamma})') ax1.plot(V_range * 1000, P_isobaric / 1e5, 'g--', linewidth=2, label=f'等圧 (P = {P1/1e5:.1f} bar)') ax1.scatter([V1 * 1000], [P1 / 1e5], color='black', s=100, zorder=5, label='初期状態') ax1.set_xlabel('Volume (L)') ax1.set_ylabel('Pressure (bar)') ax1.set_title('PV線図（可逆過程の比較）') ax1.legend() ax1.grid(True, alpha=0.3) # 仕事の計算（体積がV1→V2へ変化） V2 = 0.04 # m³ W_isothermal = n * R * T1 * np.log(V2 / V1) W_adiabatic = (P1 * V1**gamma) * (V2**(1-gamma) - V1**(1-gamma)) / (1 - gamma) W_isobaric = P1 * (V2 - V1) # 仕事の比較 ax2 = axes[1] processes = ['等温', '断熱', '等圧'] works = [W_isothermal, W_adiabatic, W_isobaric] colors = ['blue', 'red', 'green'] ax2.bar(processes, works, color=colors, alpha=0.7, edgecolor='black') ax2.set_ylabel('Work done by gas (J)') ax2.set_title(f'体積膨張の仕事（V: {V1*1000:.1f}L → {V2*1000:.1f}L）') ax2.grid(True, axis='y', alpha=0.3) for i, (process, work) in enumerate(zip(processes, works)): ax2.text(i, work + 50, f'{work:.1f} J', ha='center', fontweight='bold') plt.tight_layout() plt.savefig('thermo_reversible_processes.png', dpi=300, bbox_inches='tight') plt.show() print("=== 可逆過程における仕事 ===") print(f"初期状態: V = {V1*1000:.1f} L, P = {P1/1e5:.2f} bar, T = {T1:.1f} K") print(f"最終状態: V = {V2*1000:.1f} L\n") print(f"等温過程: W = nRT ln(V2/V1) = {W_isothermal:.2f} J") print(f"断熱過程: W = (P1V1^γ)(V2^(1-γ) - V1^(1-γ))/(1-γ) = {W_adiabatic:.2f} J") print(f"等圧過程: W = P(V2-V1) = {W_isobaric:.2f} J") print("\n→ 等温過程が最も大きな仕事を取り出せる（熱を吸収するため）")

import numpy as np import matplotlib.pyplot as plt # 物理定数 R = 8.314 # J/(mol·K) n = 1.0 # mol # 理想気体の熱力学ポテンシャル def ideal_gas_potentials(T, V, P, n, R): """理想気体の4つの熱力学ポテンシャル Args: T: 温度 (K) V: 体積 (m³) P: 圧力 (Pa) n: 物質量 (mol) R: 気体定数 (J/(mol·K)) Returns: U, H, F, G: 各ポテンシャル (J) """ # 内部エネルギー（単原子理想気体, Cv = (3/2)R） U = (3/2) * n * R * T # エンタルピー H = U + PV H = U + P * V # Helmholtz自由エネルギー F = U - TS # S = nCv ln(T) + nR ln(V) + const （Sackur-Tetrode式の簡易版） # ここでは定性的に S ∝ ln(VT^(3/2)) S = n * R * ((3/2) * np.log(T) + np.log(V) + 10) # constは適当 F = U - T * S # Gibbs自由エネルギー G = H - TS G = H - T * S return U, H, F, G, S # 温度範囲（体積固定） V_fixed = 0.0224 # m³ P_fixed = 1e5 # Pa T_range = np.linspace(100, 500, 100) potentials_vs_T = { 'U': [], 'H': [], 'F': [], 'G': [], 'S': [] } for T in T_range: U, H, F, G, S = ideal_gas_potentials(T, V_fixed, P_fixed, n, R) potentials_vs_T['U'].append(U) potentials_vs_T['H'].append(H) potentials_vs_T['F'].append(F) potentials_vs_T['G'].append(G) potentials_vs_T['S'].append(S) # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 内部エネルギー ax1 = axes[0, 0] ax1.plot(T_range, potentials_vs_T['U'], 'b-', linewidth=2) ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('U (J)') ax1.set_title('内部エネルギー U(T)') ax1.grid(True, alpha=0.3) # エンタルピー ax2 = axes[0, 1] ax2.plot(T_range, potentials_vs_T['H'], 'r-', linewidth=2) ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('H (J)') ax2.set_title('エンタルピー H(T)') ax2.grid(True, alpha=0.3) # Helmholtz自由エネルギー ax3 = axes[1, 0] ax3.plot(T_range, potentials_vs_T['F'], 'g-', linewidth=2) ax3.set_xlabel('Temperature (K)') ax3.set_ylabel('F (J)') ax3.set_title('Helmholtz自由エネルギー F(T)') ax3.grid(True, alpha=0.3) # Gibbs自由エネルギー ax4 = axes[1, 1] ax4.plot(T_range, potentials_vs_T['G'], 'purple', linewidth=2) ax4.set_xlabel('Temperature (K)') ax4.set_ylabel('G (J)') ax4.set_title('Gibbs自由エネルギー G(T)') ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_potentials.png', dpi=300, bbox_inches='tight') plt.show() # 特定温度での値 T_ref = 300 # K U, H, F, G, S = ideal_gas_potentials(T_ref, V_fixed, P_fixed, n, R) print(f"=== 理想気体の熱力学ポテンシャル（T = {T_ref} K）===") print(f"体積 V = {V_fixed * 1000:.2f} L") print(f"圧力 P = {P_fixed / 1e5:.2f} bar\n") print(f"内部エネルギー U = {U:.2f} J") print(f"エンタルピー H = U + PV = {H:.2f} J") print(f"Helmholtz自由エネルギー F = U - TS = {F:.2f} J") print(f"Gibbs自由エネルギー G = H - TS = {G:.2f} J") print(f"エントロピー S = {S:.2f} J/K") print(f"\n関係式の確認:") print(f"H - U = PV = {H - U:.2f} J (理論値: {P_fixed * V_fixed:.2f} J)") print(f"U - F = TS = {U - F:.2f} J (理論値: {T_ref * S:.2f} J)") print(f"H - G = TS = {H - G:.2f} J (理論値: {T_ref * S:.2f} J)")

import numpy as np import matplotlib.pyplot as plt # 仮想的な化学反応: A ⇌ B # ΔG = ΔH - TΔS # 反応エンタルピーとエントロピー（仮想値） delta_H = -50000 # J/mol (発熱反応) delta_S = -100 # J/(mol·K) (エントロピー減少、秩序化) def gibbs_free_energy_change(delta_H, delta_S, T): """反応のGibbs自由エネルギー変化 Args: delta_H: エンタルピー変化 (J/mol) delta_S: エントロピー変化 (J/(mol·K)) T: 温度 (K) Returns: delta_G: Gibbs自由エネルギー変化 (J/mol) """ return delta_H - T * delta_S # 温度範囲 T_range = np.linspace(200, 800, 100) delta_G_range = [gibbs_free_energy_change(delta_H, delta_S, T) for T in T_range] # 平衡温度（ΔG = 0） T_eq = delta_H / delta_S delta_G_eq = gibbs_free_energy_change(delta_H, delta_S, T_eq) # 可視化 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # ΔG vs T ax1 = axes[0] ax1.plot(T_range, delta_G_range, 'b-', linewidth=2, label='ΔG(T)') ax1.axhline(0, color='black', linestyle='--', linewidth=1.5, label='ΔG = 0 (平衡)') ax1.axvline(T_eq, color='red', linestyle='--', linewidth=1.5, label=f'平衡温度 ({T_eq:.1f} K)') ax1.fill_between(T_range, delta_G_range, 0, where=(np.array(delta_G_range) < 0), alpha=0.3, color='green', label='自発進行領域 (ΔG < 0)') ax1.fill_between(T_range, delta_G_range, 0, where=(np.array(delta_G_range) > 0), alpha=0.3, color='red', label='非自発領域 (ΔG > 0)') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('ΔG (J/mol)') ax1.set_title('Gibbs自由エネルギー変化の温度依存性') ax1.legend() ax1.grid(True, alpha=0.3) # エンタルピー項とエントロピー項の寄与 ax2 = axes[1] enthalpy_term = delta_H * np.ones_like(T_range) entropy_term = -T_range * delta_S ax2.plot(T_range, enthalpy_term, 'r-', linewidth=2, label='ΔH（エンタルピー項）') ax2.plot(T_range, entropy_term, 'b-', linewidth=2, label='-TΔS（エントロピー項）') ax2.plot(T_range, delta_G_range, 'purple', linewidth=2.5, label='ΔG = ΔH - TΔS') ax2.axhline(0, color='black', linestyle='--', linewidth=1) ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Energy (J/mol)') ax2.set_title('エンタルピー項とエントロピー項の寄与') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_gibbs_reaction.png', dpi=300, bbox_inches='tight') plt.show() print("=== 化学反応のGibbs自由エネルギー解析 ===") print(f"反応: A ⇌ B") print(f"ΔH = {delta_H / 1000:.1f} kJ/mol（発熱反応）") print(f"ΔS = {delta_S:.1f} J/(mol·K)（エントロピー減少）\n") print(f"平衡温度: T_eq = ΔH/ΔS = {T_eq:.1f} K\n") # 特定温度での判定 T_test = [300, T_eq, 700] for T in T_test: dG = gibbs_free_energy_change(delta_H, delta_S, T) if dG < 0: spontaneity = "自発的に進行（A → B）" elif dG > 0: spontaneity = "自発的に進まない（B → A が有利）" else: spontaneity = "平衡状態" print(f"T = {T:.1f} K: ΔG = {dG / 1000:.2f} kJ/mol → {spontaneity}") print("\n解釈:") print("- 低温: ΔH項が支配的 → 発熱反応が有利（ΔG < 0）") print("- 高温: -TΔS項が支配的 → エントロピー増大が重要（ΔG > 0, 逆反応有利）") print(f"- T = {T_eq:.1f} K: 両項が釣り合う → 平衡")

import numpy as np import matplotlib.pyplot as plt import sympy as sp # SymPyでのLegendre変換の例 x = sp.Symbol('x', real=True, positive=True) # 元の関数 f(x) = x² f = x**2 # df/dx df_dx = sp.diff(f, x) print("=== Legendre変換の解析例 ===") print(f"元の関数: f(x) = {f}") print(f"df/dx = {df_dx}") # p = df/dx として、x を p で表す p = sp.Symbol('p', real=True, positive=True) x_of_p = sp.solve(df_dx - p, x)[0] # x(p) を求める print(f"p = df/dx より、x(p) = {x_of_p}") # Legendre変換 g(p) = px - f(x) g = (p * x_of_p - f.subs(x, x_of_p)).simplify() print(f"Legendre変換: g(p) = px - f(x) = {g}") print() # 数値プロット x_vals = np.linspace(0.1, 3, 100) f_vals = x_vals**2 p_vals = 2 * x_vals # p = df/dx = 2x g_vals = p_vals**2 / 4 # g(p) = p²/4 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 元の関数 f(x) ax1 = axes[0] ax1.plot(x_vals, f_vals, 'b-', linewidth=2, label='f(x) = x²') ax1.set_xlabel('x') ax1.set_ylabel('f(x)') ax1.set_title('元の関数 f(x)') ax1.legend() ax1.grid(True, alpha=0.3) # Legendre変換 g(p) ax2 = axes[1] ax2.plot(p_vals, g_vals, 'r-', linewidth=2, label='g(p) = p²/4') ax2.set_xlabel('p = df/dx') ax2.set_ylabel('g(p)') ax2.set_title('Legendre変換 g(p)') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_legendre_transform.png', dpi=300, bbox_inches='tight') plt.show() # 熱力学への応用例 print("=== 熱力学ポテンシャルへの応用 ===") print("内部エネルギー U(S,V):") print(" dU = TdS - PdV") print(" → T = (∂U/∂S)_V, P = -(∂U/∂V)_S") print() print("Legendre変換:") print(" S → T: F = U - TS （Helmholtz自由エネルギー）") print(" V → P: H = U + PV （エンタルピー）") print(" 両方: G = U - TS + PV （Gibbs自由エネルギー）") print() print("この変換により、実験で制御しやすい変数（T, P）を") print("自然変数とするポテンシャルを構築できる。")

import numpy as np import matplotlib.pyplot as plt # Debye理論によるエントロピーと比熱 def debye_entropy(T, theta_D, R): """Debye理論によるエントロピー（簡易版） Args: T: 温度 (K) theta_D: Debye温度 (K) R: 気体定数 (J/(mol·K)) Returns: S: エントロピー (J/(mol·K)) """ x = theta_D / T if T < 0.01: # 数値計算の安定性 return 0 # 低温極限: S ∝ (T/θ_D)³ if x > 10: return (12/5) * np.pi**4 * R * (T / theta_D)**3 # 高温極限: S → 3R ln(T/θ_D) + const else: return 3 * R * (np.log(T / theta_D) + 4/3) def debye_heat_capacity(T, theta_D, R): """Debye理論による定積比熱 Args: T: 温度 (K) theta_D: Debye温度 (K) R: 気体定数 (J/(mol·K)) Returns: C_V: 定積比熱 (J/(mol·K)) """ x = theta_D / T # 低温極限: C_V ∝ T³ if x > 10: return (12/5) * np.pi**4 * R * (T / theta_D)**3 # 高温極限: C_V → 3R (Dulong-Petit) else: return 3 * R # 物理定数 R = 8.314 # J/(mol·K) theta_D_Cu = 343 # Cuのebye温度 (K) # 温度範囲 T_range = np.logspace(-1, 3, 200) # 0.1 K to 1000 K S_vals = [debye_entropy(T, theta_D_Cu, R) for T in T_range] C_vals = [debye_heat_capacity(T, theta_D_Cu, R) for T in T_range] # 可視化 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # エントロピー ax1 = axes[0] ax1.loglog(T_range, S_vals, 'b-', linewidth=2) ax1.axvline(theta_D_Cu, color='red', linestyle='--', linewidth=1.5, label=f'Debye温度 ({theta_D_Cu} K)') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Entropy (J/(mol·K))') ax1.set_title('エントロピーの温度依存性（Cu）') ax1.legend() ax1.grid(True, alpha=0.3, which='both') # 比熱 ax2 = axes[1] ax2.loglog(T_range, C_vals, 'r-', linewidth=2, label='Debye理論') ax2.axhline(3 * R, color='green', linestyle='--', linewidth=1.5, label='Dulong-Petit則 (3R)') ax2.axvline(theta_D_Cu, color='black', linestyle='--', linewidth=1.5, label=f'Debye温度') ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Heat capacity C_V (J/(mol·K))') ax2.set_title('定積比熱の温度依存性（Cu）') ax2.legend() ax2.grid(True, alpha=0.3, which='both') plt.tight_layout() plt.savefig('thermo_third_law.png', dpi=300, bbox_inches='tight') plt.show() print("=== 熱力学第3法則の検証（Cu）===") print(f"Debye温度: θ_D = {theta_D_Cu} K\n") # 低温でのエントロピー T_low = [0.1, 1, 10, 100] print("低温でのエントロピー:") for T in T_low: S = debye_entropy(T, theta_D_Cu, R) print(f" T = {T:6.1f} K: S = {S:.6f} J/(mol·K)") print(f"\n→ T → 0 で S → 0 （第3法則）") print(f"\n低温での振る舞い:") print(f" S ∝ T³ (T << θ_D)") print(f" C_V ∝ T³ (T << θ_D)") print(f"\n高温での振る舞い:") print(f" C_V → 3R = {3*R:.2f} J/(mol·K) (Dulong-Petit則)")

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve # van der Waals定数（N₂） a_N2 = 0.1408 # Pa·m⁶/mol² b_N2 = 3.913e-5 # m³/mol R = 8.314 # J/(mol·K) C_V = 2.5 * R # J/(mol·K)（二原子分子） def van_der_waals_pressure(V, n, T, a, b, R): """van der Waals圧力""" return n * R * T / (V - n * b) - a * n**2 / V**2 def van_der_waals_internal_energy(V, n, T, a, C_V): """van der Waals内部エネルギー""" return n * C_V * T - a * n**2 / V # 固定パラメータ n = 1.0 # mol T = 300 # K # 体積範囲 V_range = np.linspace(0.001, 0.1, 200) # van der Waals圧力と内部エネルギー P_vdw = [van_der_waals_pressure(V, n, T, a_N2, b_N2, R) for V in V_range] U_vdw = [van_der_waals_internal_energy(V, n, T, a_N2, C_V) for V in V_range] # 理想気体との比較 P_ideal = [n * R * T / V for V in V_range] U_ideal = n * C_V * T * np.ones_like(V_range) # 可視化 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 圧力比較 ax1 = axes[0] ax1.plot(V_range * 1000, np.array(P_vdw) / 1e5, 'b-', linewidth=2, label='van der Waals') ax1.plot(V_range * 1000, np.array(P_ideal) / 1e5, 'r--', linewidth=2, label='理想気体') ax1.set_xlabel('Volume (L)') ax1.set_ylabel('Pressure (bar)') ax1.set_title(f'PV等温線（T = {T} K, N₂）') ax1.set_xlim([0, 100]) ax1.set_ylim([0, 50]) ax1.legend() ax1.grid(True, alpha=0.3) # 内部エネルギー比較 ax2 = axes[1] ax2.plot(V_range * 1000, U_vdw, 'b-', linewidth=2, label='van der Waals') ax2.plot(V_range * 1000, U_ideal, 'r--', linewidth=2, label='理想気体') ax2.set_xlabel('Volume (L)') ax2.set_ylabel('Internal energy U (J)') ax2.set_title(f'内部エネルギー（T = {T} K, N₂）') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_van_der_waals.png', dpi=300, bbox_inches='tight') plt.show() # 数値比較 V_test = 0.0224 # m³ (標準状態) P_vdw_val = van_der_waals_pressure(V_test, n, T, a_N2, b_N2, R) P_ideal_val = n * R * T / V_test U_vdw_val = van_der_waals_internal_energy(V_test, n, T, a_N2, C_V) U_ideal_val = n * C_V * T print("=== van der Waals気体 vs 理想気体（N₂, T = 300 K）===") print(f"体積 V = {V_test * 1000:.2f} L\n") print("圧力:") print(f" van der Waals: P = {P_vdw_val / 1e5:.4f} bar") print(f" 理想気体: P = {P_ideal_val / 1e5:.4f} bar") print(f" 差: {(P_vdw_val - P_ideal_val) / 1e5:.4f} bar\n") print("内部エネルギー:") print(f" van der Waals: U = {U_vdw_val:.2f} J") print(f" 理想気体: U = {U_ideal_val:.2f} J") print(f" 差: {U_vdw_val - U_ideal_val:.2f} J") print(f"\n補正項:") print(f" 引力補正: -an²/V = {-a_N2 * n**2 / V_test:.2f} J") print(f" 体積補正: nb = {n * b_N2 * 1000:.4f} L")

import numpy as np import matplotlib.pyplot as plt # Grüneisen関係式: α = γC_V / (VK_T) # α: 熱膨張率, γ: Grüneisenパラメータ, K_T: 等温圧縮率 def thermal_expansion_coefficient(T, gamma, C_V, V, K_T): """熱膨張率 α = γC_V/(VK_T)""" return gamma * C_V(T) / (V * K_T) def debye_heat_capacity_solid(T, theta_D, R): """固体のDebye比熱""" x = theta_D / T if x > 10: # 低温 return (12/5) * np.pi**4 * R * (T / theta_D)**3 else: # 高温（Dulong-Petit） return 3 * R # Cuの物性値 theta_D_Cu = 343 # K gamma_Cu = 2.0 # Grüneisenパラメータ V_Cu = 7.11e-6 # m³/mol K_T_Cu = 1.4e11 # Pa R = 8.314 # J/(mol·K) # 温度範囲 T_range = np.linspace(10, 1000, 200) # 熱膨張率 alpha_vals = [] C_V_vals = [] for T in T_range: C_V = debye_heat_capacity_solid(T, theta_D_Cu, R) alpha = thermal_expansion_coefficient(T, gamma_Cu, lambda t: C_V, V_Cu, K_T_Cu) alpha_vals.append(alpha) C_V_vals.append(C_V) # 可視化 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 熱膨張率 ax1 = axes[0] ax1.plot(T_range, np.array(alpha_vals) * 1e6, 'b-', linewidth=2) ax1.axvline(theta_D_Cu, color='red', linestyle='--', linewidth=1.5, label=f'Debye温度 ({theta_D_Cu} K)') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Thermal expansion coefficient α (10⁻⁶ K⁻¹)') ax1.set_title('熱膨張率の温度依存性（Cu）') ax1.legend() ax1.grid(True, alpha=0.3) # 比熱との相関 ax2 = axes[1] ax2.scatter(C_V_vals, np.array(alpha_vals) * 1e6, c=T_range, cmap='viridis', s=20) cbar = plt.colorbar(ax2.scatter(C_V_vals, np.array(alpha_vals) * 1e6, c=T_range, cmap='viridis', s=20), ax=ax2) cbar.set_label('Temperature (K)') ax2.set_xlabel('Heat capacity C_V (J/(mol·K))') ax2.set_ylabel('Thermal expansion coefficient α (10⁻⁶ K⁻¹)') ax2.set_title('熱膨張率と比熱の関係（Grüneisen関係式）') ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_thermal_expansion.png', dpi=300, bbox_inches='tight') plt.show() print("=== 固体の熱膨張（Cu）===") print(f"Debye温度: θ_D = {theta_D_Cu} K") print(f"Grüneisenパラメータ: γ = {gamma_Cu}") print(f"等温圧縮率: K_T = {K_T_Cu:.2e} Pa\n") # 特定温度での値 T_test = [100, 300, 500, 1000] print("温度依存性:") for T in T_test: C_V = debye_heat_capacity_solid(T, theta_D_Cu, R) alpha = thermal_expansion_coefficient(T, gamma_Cu, lambda t: C_V, V_Cu, K_T_Cu) print(f"T = {T:4d} K: C_V = {C_V:.2f} J/(mol·K), α = {alpha*1e6:.2f} × 10⁻⁶ K⁻¹") print("\nGrüneisen関係式: α = γC_V/(VK_T)") print("→ 熱膨張率は比熱に比例（低温で両方ともゼロに近づく）") print("→ 材料設計における熱応力評価に重要")