問題
化学反応シミュレーションや燃焼計算の分野では、各化学種の熱力学関数(定圧モル熱容量 \(C_p\)、エンタルピー \(H\)、エントロピー \(S\))を温度の関数として表すために、NASA多項式形式(7係数多項式)が用いられることがあります。
この形式では、無次元化された熱力学関数が次のように表されます。
\[\frac{C_p}{R} = a_1 + a_2 T + a_3 T^2 + a_4 T^3 + a_5 T^4\]
\[\frac{H}{RT} = a_1 + \frac{a_2}{2}T + \frac{a_3}{3}T^2 + \frac{a_4}{4}T^3 + \frac{a_5}{5}T^4 + \frac{a_6}{T}\]
\[\frac{S}{R} = a_1 \ln T + a_2 T + \frac{a_3}{2}T^2 + \frac{a_4}{3}T^3 + \frac{a_5}{4}T^4 + a_7\]
ここで \(R\) は気体定数、\(a_1,\ldots,a_7\) は各化学種・各温度範囲ごとに定められた無次元の係数です。
\(H/RT\) の式は \(C_p/R\) を温度で積分することで(積分定数として \(a_6/T\) の項が現れる形で)導かれ、\(S/R\) の式は \(C_p/(RT)\) を積分することで導かれる、という関係でつながっています。
窒素分子 N₂(g) の 300〜1000 K の範囲における係数として、次の値が与えられているとします。
| 係数 | 値 |
|---|---|
| \(a_1\) | \(3.53100528\) |
| \(a_2\) | \(-1.23660988\times10^{-4}\) |
| \(a_3\) | \(-5.02999433\times10^{-7}\) |
| \(a_4\) | \(2.43530612\times10^{-9}\) |
| \(a_5\) | \(-1.40881235\times10^{-12}\) |
| \(a_6\) | \(-1.04697628\times10^{3}\) |
| \(a_7\) | \(2.96747038\) |
(1) これらの係数を用いて、\(T=300\sim1000~\text{K}\) における \(C_p(T)\)、\(H(T)\)、\(S(T)\) をRで計算し、可視化しなさい。
(2) \(C_p/R\) の式を \(H/RT\) の式から数値微分(\(\frac{dH}{dT}=C_p\) の関係を利用)によって逆算し、多項式から直接計算した \(C_p\) と一致することを確認しなさい。
(3) \(T=298.15~\text{K}\) における \(S\) の値を求め、これが窒素分子の標準モルエントロピー(文献値 約191.6 J/(mol·K))と近い値になることを確認しなさい。
解答方針
NASA多項式は「\(C_p/R\)、\(H/RT\)、\(S/R\)」という3つの無次元関数を、共通の係数 \(a_1,\ldots,a_7\) から生成する仕組みになっています。
Rコードで、これら3つの式をそれぞれ関数として実装し、共通の係数ベクトルを渡すことで一括して計算できるようにします。
(2)の検証では、\(H(T)\) を数値的に温度微分した結果が、多項式から直接得られる \(C_p(T)\) と一致するかを確認することで、NASA多項式の内部的な整合性(\(C_p\)と\(H\)が同じ係数から矛盾なく導かれていること)を検証します。
# ------------------------------------------------------------
# NASA多項式形式による熱力学関数の計算: N2(g), 300-1000K
# ------------------------------------------------------------
R_gas <- 8.314 # 気体定数 J/(mol・K)
# N2(g)のNASA多項式係数(300-1000K区間)
a_coef <- c(
a1 = 3.53100528,
a2 = -1.23660988e-4,
a3 = -5.02999433e-7,
a4 = 2.43530612e-9,
a5 = -1.40881235e-12,
a6 = -1.04697628e3,
a7 = 2.96747038
)
# Cp/R を計算する関数
Cp_over_R <- function(T, a) {
a["a1"] + a["a2"] * T + a["a3"] * T^2 + a["a4"] * T^3 + a["a5"] * T^4
}
# H/RT を計算する関数
H_over_RT <- function(T, a) {
a["a1"] + (a["a2"] / 2) * T + (a["a3"] / 3) * T^2 +
(a["a4"] / 4) * T^3 + (a["a5"] / 5) * T^4 + a["a6"] / T
}
# S/R を計算する関数
S_over_R <- function(T, a) {
a["a1"] * log(T) + a["a2"] * T + (a["a3"] / 2) * T^2 +
(a["a4"] / 3) * T^3 + (a["a5"] / 4) * T^4 + a["a7"]
}
# 次元付きの熱力学関数として計算する関数群
Cp_func <- function(T, a) R_gas * Cp_over_R(T, a)
H_func <- function(T, a) R_gas * T * H_over_RT(T, a)
S_func <- function(T, a) R_gas * S_over_R(T, a)
# --- (1) 300~1000Kにおける計算と可視化 ---
T_seq <- seq(300, 1000, by = 10)
Cp_vals <- sapply(T_seq, Cp_func, a = a_coef)
H_vals <- sapply(T_seq, H_func, a = a_coef) / 1000 # kJ/molに変換
S_vals <- sapply(T_seq, S_func, a = a_coef)
df_nasa <- data.frame(
T_K = T_seq,
Cp_J_molK = Cp_vals,
H_kJ_mol = H_vals,
S_J_molK = S_vals
)
print(head(df_nasa, 10)) T_K Cp_J_molK H_kJ_mol S_J_molK
1 300 29.12377 0.05387785 191.7781
2 310 29.13119 0.34515133 192.7332
3 320 29.14019 0.63650689 193.6582
4 330 29.15080 0.92796050 194.5551
5 340 29.16306 1.21952841 195.4255
6 350 29.17698 1.51122718 196.2710
7 360 29.19260 1.80307364 197.0932
8 370 29.20993 2.09508482 197.8933
9 380 29.22899 2.38727797 198.6725
10 390 29.24981 2.67967051 199.4320可視化
library(ggplot2)
library(patchwork) # 複数プロットを並べて表示するために使用
p_cp <- ggplot(df_nasa, aes(x = T_K, y = Cp_J_molK)) +
geom_line(color = "darkorange", linewidth = 1) +
labs(title = expression(C[p](T)), x = "T (K)", y = "Cp (J/(mol・K))") +
theme_minimal(base_size = 11)
p_h <- ggplot(df_nasa, aes(x = T_K, y = H_kJ_mol)) +
geom_line(color = "steelblue", linewidth = 1) +
labs(title = "H(T)", x = "T (K)", y = "H (kJ/mol)") +
theme_minimal(base_size = 11)
p_s <- ggplot(df_nasa, aes(x = T_K, y = S_J_molK)) +
geom_line(color = "seagreen", linewidth = 1) +
labs(title = "S(T)", x = "T (K)", y = "S (J/(mol・K))") +
theme_minimal(base_size = 11)
p_combined <- p_cp / p_h / p_s # 縦に3段構成で表示
print(p_combined)
(2) H(T)の数値微分によるCp(T)の逆算と検証
# H(T)を数値的に温度微分してCp(T)を逆算する
# 中心差分法: dH/dT ≈ (H(T+h) - H(T-h)) / (2h)
numeric_derivative <- function(f, T, a, h = 0.01) {
(f(T + h, a) - f(T - h, a)) / (2 * h)
}
Cp_from_H_numeric <- sapply(T_seq, function(Ti) {
numeric_derivative(H_func, Ti, a_coef)
})
Cp_direct <- sapply(T_seq, Cp_func, a = a_coef)
compare_df <- data.frame(
T_K = T_seq,
Cp_direct = round(Cp_direct, 6),
Cp_from_H_deriv = round(Cp_from_H_numeric, 6),
diff = round(Cp_direct - Cp_from_H_numeric, 8)
)
cat("=== Cp(T)直接計算 vs H(T)の数値微分による逆算 の比較(一部抜粋)===\n")
print(compare_df[seq(1, nrow(compare_df), by = 20), ])=== Cp(T)直接計算 vs H(T)の数値微分による逆算 の比較(一部抜粋)===
T_K Cp_direct Cp_from_H_deriv diff
1 300 29.12377 29.12377 0
21 500 29.59607 29.59607 0
41 700 30.72046 30.72046 0
61 900 32.11945 32.11945 0数値微分による逆算値と多項式からの直接計算値は(数値誤差の範囲で)完全に一致しており、\(H/RT\) の式が確かに \(C_p/R\) を積分した形になっていることが確認できます。
(3) 298.15 K における標準モルエントロピーとの比較
T_std <- 298.15
S_298 <- S_func(T_std, a_coef)
cat("=== 298.15 K における N2(g) のエントロピー ===\n")
cat(sprintf("NASA多項式からの計算値: S = %.2f J/(mol・K)\n", S_298))
cat("文献の標準モルエントロピー(参考値): 約191.6 J/(mol・K)\n")
cat(sprintf("差: %.2f J/(mol・K)\n", S_298 - 191.6))=== 298.15 K における N2(g) のエントロピー ===
NASA多項式からの計算値: S = 191.60 J/(mol・K)
文献の標準モルエントロピー(参考値): 約191.6 J/(mol・K)
差: -0.00 J/(mol・K)計算値は文献値とほぼ一致しており、NASA多項式が実験データに基づいて精度よくフィッティングされていることが確認できます。
まとめ
| 項目 | 内容 |
|---|---|
| NASA多項式形式 | \(C_p/R,\ H/RT,\ S/R\) を共通の7係数(\(a_1\sim a_7\))で表現する形式 |
| 係数間の関係 | \(H/RT\) の式は \(C_p/R\) を積分した形、\(S/R\) の式は \(C_p/(RT)\) を積分した形になっている |
| 検証(2) | \(H(T)\) の数値微分から \(C_p(T)\) を逆算し、直接計算値と一致することを確認 |
| 検証(3) | 298.15 K での計算値が文献の標準モルエントロピーとよく一致 |
補足
標準モルエントロピーとは
標準モルエントロピーとは、ある物質1molが標準状態(通常、圧力1 bar、指定された温度、最も安定な状態)にあるときのエントロピーの絶対値のことで、記号 \(S^\circ\) で表されます(単位 J/(mol·K))。
エントロピー \(S\) は系の乱雑さ・とりうる微視的状態の多さを表す状態量であり、統計力学的には
\[S = k_B \ln W\]
という関係(ボルツマンの関係式)で微視的な描像と結びついています。ここで \(k_B\) はボルツマン定数(単位 J/K)、\(W\) はその巨視的状態に対応する微視的状態の数(無次元の場合の数)です。
熱力学第三法則により、完全結晶の絶対零度(\(T=0~\text{K}\))におけるエントロピーは0と定義されるため、そこを基準として温度を上げていったときのエントロピー変化を積算すれば、任意の温度における絶対的なエントロピー値 \(S^\circ\) を決定できます。
この値は、反応エントロピー \(\Delta S^\circ\) を生成物と反応物の \(S^\circ\) の差として計算する際の基礎データとして用いられます。
なぜモルで表すのか
物質は原子・分子という極めて小さな粒子からできており、実験室で扱う量(グラムやリットル単位)の中には天文学的な数の粒子が含まれています。例えば水18gの中には水分子が
\[N = n \times N_A\]
個含まれます。ここで \(n\) は物質量(単位 mol)、\(N_A\) はアボガドロ定数(1molあたりの粒子数、単位 mol⁻¹)、\(N\) は粒子の個数です。
このように桁数の大きい \(N\) を直接扱う代わりに、\(n\)という扱いやすい数値に置き換えられる点が、モルを使う最大の利点です。
さらに、化学反応式の係数はそのまま物質量の比を表しているため、\(n_A/n_B\) のようにモル単位で量を扱えば、反応に関わる物質同士の量的関係(化学量論)を係数の比だけで直接計算できます。
質量や体積は物質ごとに換算が必要ですが、モルを介せば異なる物質同士の量を統一的に比較・計算できるという点が、実験や理論計算の両面で大きな利便性をもたらしています。
質量から物質量への変換方法
物質の質量(グラム)から物質量(モル)へ変換するには、モル質量を用いた次の関係式を使います。
\[n = \frac{m}{M}\]
ここで \(n\) は物質量(単位 mol)、\(m\) は物質の質量(単位 g)、\(M\) はモル質量(その物質1molあたりの質量、単位 g/mol)です。
モル質量 \(M\) は、物質を構成する原子の原子量(周期表に示されている、各元素1molあたりのおおよその質量に対応する数値)を、化学式に含まれる原子の数だけ足し合わせることで求められます。
例えば水(\(\text{H}_2\text{O}\))であれば、水素の原子量(約1.0)を2個分、酸素の原子量(約16.0)を1個分足し合わせて
\[M_{\text{H}_2\text{O}} = 2\times1.0 + 1\times16.0 = 18.0~\text{g/mol}\]
となります。
したがって、例えば水36gの物質量を求めたい場合は、\(n = 36/18.0 = 2.0~\text{mol}\) という計算で求められます。
この変換により、実験室で測定できる質量というマクロな量から、粒子の個数という物質量の指標へと橋渡しすることができます。
完全結晶とは
完全結晶とは、原子や分子が結晶格子の中で一切の乱れなく、完全に規則正しい周期構造で配列した理想的な結晶のことを指します。
格子欠陥・不純物の混入・分子の配向の乱れなどが一切存在しない、理論上の理想状態です。
熱力学第三法則では、絶対零度(\(T=0~\text{K}\))における完全結晶のエントロピーは0であると定められています。これは
\[S = k_B \ln W\]
という関係式(ここで \(S\) はエントロピー、単位 J/K。\(k_B\) はボルツマン定数、単位 J/K。\(W\) はその状態を実現する微視的配置の数、無次元の場合の数)において、完全結晶かつ絶対零度の状態では原子・分子の配置が一意に定まり、\(W=1\) となる(\(\ln 1 = 0\))ことに対応しています。
この基準があることで、任意の温度・物質における絶対的なエントロピー値 \(S^\circ\)(標準モルエントロピー)を、絶対零度からの変化量を積算する形で一意に定義・測定できるようになります。
実際の結晶には格子欠陥などの乱れが存在するため、完全結晶はあくまで理論上の基準として位置づけられます。
以上です。
