Rで理論化学:キルヒホッフの法則(反応エンタルピーの温度依存性)

問題

アンモニア合成反応(ハーバー・ボッシュ法)

\[\text{N}_2(g) + 3\text{H}_2(g) \rightleftharpoons 2\text{NH}_3(g)\]

の標準反応エンタルピーは 298 K において

\[\Delta H^\circ_{298} = -92.4~\text{kJ/mol}\]

であるとします。

各成分の定圧モル熱容量 \(C_p\)(単位: J/(mol·K))は、次の温度依存性を持つ経験式(温度 \(T\) は K 単位)で近似できるとします。

(1) キルヒホッフの法則

\[\left(\frac{\partial \Delta H}{\partial T}\right)_p = \Delta C_p\]

に基づき、\(\Delta C_p(T) = 2C_{p,\text{NH}_3} - C_{p,\text{N}_2} - 3C_{p,\text{H}_2}\) をRで多項式として構成し、\(\Delta H(T)\) を \(298~\text{K}\) を基準として

\[\Delta H(T) = \Delta H^\circ_{298} + \int_{298}^{T} \Delta C_p(T')\, dT'\]

により、\(T = 200 \sim 800~\text{K}\) の範囲で計算しなさい。

積分はRの数値積分(integrate())を用いる方法と、多項式を解析的に積分する方法の両方で行い、結果が一致することを確認しなさい。

(2) 計算結果を可視化し、温度が上がるにつれて \(\Delta H\) がどのように変化するか(発熱がより強くなる/弱くなる)を考察しなさい。

解答方針

まず各物質の \(C_p(T)\) を係数ベクトルとしてRに実装し、反応の \(\Delta C_p(T)\) を線形結合として求めます。

次に、\(\Delta H(T)\) を「数値積分」と「解析的な多項式積分」の2通りで計算して整合性を検証します。

最後に、得られた \(\Delta H(T)\) を温度に対してプロットします。

# ------------------------------------------------------------
# キルヒホッフの法則:反応エンタルピーの温度依存性
# N2(g) + 3H2(g) -> 2NH3(g)
# ------------------------------------------------------------

# 各物質のCp(T) = a + b*T + c*T^2 の係数(J/(mol・K))
Cp_coef <- list(
  N2  = c(a = 27.9, b = 4.27e-3, c = 0),
  H2  = c(a = 27.3, b = 3.30e-3, c = 0.50e-6),
  NH3 = c(a = 25.9, b = 33.0e-3, c = -3.05e-6)
)

# Cp(T)を計算する汎用関数
Cp_func <- function(T, coef) {
  coef["a"] + coef["b"] * T + coef["c"] * T^2
}

# 反応の化学量論係数: 2 NH3 - N2 - 3 H2
nu <- c(N2 = -1, H2 = -3, NH3 = 2)

# 反応のΔCp(T)を計算する関数(各成分のCp(T)を化学量論係数で線形結合)
dCp_func <- function(T) {
  nu["N2"] * Cp_func(T, Cp_coef$N2) +
    nu["H2"] * Cp_func(T, Cp_coef$H2) +
    nu["NH3"] * Cp_func(T, Cp_coef$NH3)
}

# ΔCp(T)の多項式係数を直接合成しておく(解析的積分に使う)
dCp_a <- nu["N2"] * Cp_coef$N2["a"] + nu["H2"] * Cp_coef$H2["a"] + nu["NH3"] * Cp_coef$NH3["a"]
dCp_b <- nu["N2"] * Cp_coef$N2["b"] + nu["H2"] * Cp_coef$H2["b"] + nu["NH3"] * Cp_coef$NH3["b"]
dCp_c <- nu["N2"] * Cp_coef$N2["c"] + nu["H2"] * Cp_coef$H2["c"] + nu["NH3"] * Cp_coef$NH3["c"]

cat("=== ΔCp(T) の多項式係数 ===\n")
cat(sprintf("a = %.4f, b = %.6e, c = %.6e\n\n", dCp_a, dCp_b, dCp_c))

# 298 KにおけるΔCpを確認(参考値)
cat(sprintf("ΔCp(298 K) = %.3f J/(mol・K)\n\n", dCp_func(298)))

=== ΔCp(T) の多項式係数 ===
a = -58.0000, b = 5.183000e-02, c = -7.600000e-06

ΔCp(298 K) = -43.230 J/(mol・K)

dH_298 <- -92.4e3 # J/mol (kJ/molをJ/molに変換)
T0 <- 298

# --- 方法(a): Rの数値積分 integrate() を使う ---
dH_numeric <- function(T) {
  sapply(T, function(Ti) {
    integral <- integrate(dCp_func, lower = T0, upper = Ti)$value
    dH_298 + integral
  })
}

# --- 方法(b): 多項式の解析的積分 ---
# ∫(a + bT + cT^2)dT = a*T + b/2*T^2 + c/3*T^3 + const
# 定積分: [a*T + b/2*T^2 + c/3*T^3] を T0 から T まで評価
antideriv <- function(T) {
  dCp_a * T + (dCp_b / 2) * T^2 + (dCp_c / 3) * T^3
}

dH_analytic <- function(T) {
  dH_298 + (antideriv(T) - antideriv(T0))
}

# 200K~800Kの範囲で両者を比較
T_seq <- seq(200, 800, by = 50)

dH_num <- dH_numeric(T_seq)
dH_ana <- dH_analytic(T_seq)

comparison <- data.frame(
  T_K = T_seq,
  dH_numeric_kJ = round(dH_num / 1000, 4),
  dH_analytic_kJ = round(dH_ana / 1000, 4),
  diff = round((dH_num - dH_ana), 8)
)

print(comparison)

   T_K dH_numeric_kJ dH_analytic_kJ diff
1  200      -87.9340       -87.9340    0
2  250      -90.2702       -90.2702    0
3  300      -92.4864       -92.4864    0
4  350      -94.5843       -94.5843    0
5  400      -96.5660       -96.5660    0
6  450      -98.4334       -98.4334    0
7  500     -100.1882      -100.1882    0
8  550     -101.8325      -101.8325    0
9  600     -103.3681      -103.3681    0
10 650     -104.7969      -104.7969    0
11 700     -106.1209      -106.1209    0
12 750     -107.3419      -107.3419    0
13 800     -108.4618      -108.4618    0

可視化

library(ggplot2)

T_fine <- seq(200, 800, length.out = 300)
dH_fine <- dH_analytic(T_fine) / 1000 # kJ/molに変換

df_plot <- data.frame(T_K = T_fine, dH_kJ = dH_fine)

p <- ggplot(df_plot, aes(x = T_K, y = dH_kJ)) +
  geom_line(color = "firebrick", linewidth = 1) +
  geom_point(
    data = data.frame(T_K = 298, dH_kJ = -92.4),
    color = "black", size = 3
  ) +
  annotate("text",
    x = 298, y = -92.4, label = "基準点 (298 K, -92.4 kJ/mol)",
    hjust = -0.1, vjust = -0.8, size = 3.5
  ) +
  labs(
    title = "アンモニア合成反応のΔHの温度依存性(キルヒホッフの法則)",
    x = "温度 T (K)",
    y = expression(Delta * H ~ "(kJ/mol)")
  ) +
  theme_minimal(base_size = 13)

print(p)

---

考察

計算結果から、\(\Delta C_p(298~\text{K}) \approx -43.2~\text{J/(mol·K)}\) と、負でかつ比較的大きな絶対値を持つことが分かります。

これは、生成物側(NH₃ 2 mol)の熱容量の合計が、反応物側(N₂ 1 mol + H₂ 3 mol)の熱容量の合計よりもかなり小さいことを意味しています。

\(\Delta C_p < 0\) であることから、キルヒホッフの法則 \(\left(\dfrac{\partial \Delta H}{\partial T}\right)_p = \Delta C_p\) より、温度が上昇するにつれて \(\Delta H\) は単調に減少し続けます(発熱がより強くなる方向へ変化します)。

実際、計算結果でも \(T=200~\text{K}\) で約 \(-87.9~\text{kJ/mol}\) であった \(\Delta H\) が、\(T=800~\text{K}\) では約 \(-108.5~\text{kJ/mol}\) まで減少しており、200 Kから800 Kの範囲でおよそ20.5 kJ/molもの変化が生じています。

この変化幅は、基準となる \(\Delta H^\circ_{298}=-92.4~\text{kJ/mol}\) の絶対値と比較して約22%に相当し、決して無視できる大きさではありません。

したがって、この反応系においては\(\Delta H\)を温度に依存しない定数として扱う近似は適切ではなく、特に高温条件(工業的なアンモニア合成は400〜500℃程度の高温で行われることが多い)での平衡計算や熱収支の見積もりを行う際には、キルヒホッフの法則に基づいて \(\Delta C_p(T)\) の温度依存性を明示的に考慮する必要があることが、この結果から示唆されます。

まとめ

補足