メタン水蒸気改質反応の平衡計算

出村雅彦

2011年12月22日

取り扱う問題とアプローチ

反応管にメタンを25 ml/min、水蒸気を25 ml/min流そう。反応管を加熱するとメタン水蒸気改質反応や水性シフト反応が起こる。化学的平衡に達すると、出口での各ガスの流量、そして、メタン転化率はどうなるか？--これが、ここで取り扱う問題である。メタン・水蒸気比（S/C）や不活性バッファーガス（N

等）導入が、平衡メタン転化率にどんな影響を及ぼすのかも、視野のうちに入れたい。

平衡状態では、それ以上反応が進行せず、Gibbsエネルギーが極小化する。Gibbsエネルギーが極小化する反応の進行度合いを求めれば良いわけだ。その方法には、二つある。一つは、平衡定数から求める方法。これは、昔からある古典的な方法である。ここでは、この方法について、詳しく述べる。もう一つの方法は、Gibbsエネルギーを直接最小化する方法である。これは平衡の定義そのものである。この方法は、フリーで公開されているパッケージ（CEA）を利用すると簡単に試すことができる。ここでは、その使い方を紹介するだけにとどめ、詳しい解説はしない（というか、勉強不足でできない）。もちろん、基本的にどちらの方法も、同じ答えを与える。

では、まずは、平衡定数から計算する方法から、見てみよう。古典な方法ではあるが、これを理解しておくと、計算結果をより深く理解できると思う。

平衡定数から計算する方法

平衡定数と平衡ガス分圧

ある反応の平衡定数は、反応式右辺（反応物）のガスの平衡分圧を分母に、左辺（生成物）の平衡ガス分圧を分母においた形をしている。メタン水蒸気改質反応で説明しよう。反応式は、以下で表される。

$\begin{displaymath} \ensuremath{\xi_{\mathrm{M}}}: \ensuremath{\mathrm{CH}_4}+ \... ...thrm{O}}= \ensuremath{\mathrm{CO}}+ 3\ensuremath{\mathrm{H}_2} \end{displaymath}$

(1)

$\ensuremath{\xi_{\mathrm{M}}}$ は反応の進行度を表す変数で、後ほど登場する。この反応式で書いた時のメタン水蒸気改質反応の平衡定数は、

$\displaystyle \ensuremath{K_{\mathrm{M}}}$	$\textstyle =$	$\displaystyle \frac{\frac{P_{\ensuremath{\mathrm{CO}}}}{P^0} \cdot \left(\frac{... ...mathrm{CH}_4}}}{P^0} \cdot \frac{P_{\ensuremath{\mathrm{H}_2\mathrm{O}}}}{P^0}}$	(2)
$\displaystyle \ensuremath{K_{\mathrm{M}}}$	$\textstyle =$	$\displaystyle \frac{P_{\ensuremath{\mathrm{CO}}}\cdot{P_{\ensuremath{\mathrm{H}... ...}{P_{\ensuremath{\mathrm{CH}_4}} \cdot P_{\ensuremath{\mathrm{H}_2\mathrm{O}}}}$	(3)

となる。ここで、

（

は物質名）は、平衡に達した時の部分モル分圧（以下、ガス分圧と呼ぶ）、

は標準状態の圧力で通常は1 barである。二つ式があるが、

barなので、計算する上で違いはない。以下では、式（3）を採用する。ただし、各ガスのガス分圧は1 barに対する比で表されていると考えて、 $\ensuremath{K_{\mathrm{M}}}$ は無次元数であることを忘れないようにしたい。

水素のガス分圧が3乗となっていることに注意してほしい。反応式中の係数でベキ乗するわけである。ということは、式（1）を2倍すると、その平衡定数は式（3）を2乗したものになる。つまり、反応式を与えないと平衡定数は決まらない。メタン水蒸気改質反応の場合はわざわざ2倍した式を使うことはないが、一酸化炭素の酸化反応（COシフト反応）などでは、一酸化炭素の係数を1にして書いてみたり、分数が出てくるのを避けるために係数を2にしてみたり、目的によって反応式も変わり、それに伴って、平衡定数の値も変わることになるので、注意が必要である。

平衡状態を求めるということ

平衡定数は、（純粋な）各ガス分子の標準状態（1 bar）のGibbs自由エネルギー（化学ポテンシャル）から求めることができ、反応と温度を決めれば、一意にその値が決まる。その値に合うように、各ガス分子のガス分圧を決めれば良いわけである。これが決まれば、平衡状態の流量が計算でき、メタンの平衡転化率や各ガスの平衡選択率が計算できる。

例えば、800のメタン水蒸気改質反応の平衡定数は、 $\ensuremath{K_{\mathrm{M}}}\simeq 1.57 \cdot 10^2$ である。この値になるように、ガス分圧を決めてみよう。 $P_{\ensuremath{\mathrm{CO}}} = 1.57 \cdot 10^2$ barで、他が全て1 barというガス分圧の組み合わせは、平衡定数と一致する。従って、ひとつの平衡状態である。似たような解として、 $P_{\ensuremath{\mathrm{H}_2}}=\sqrt[3]{1.57 \cdot 10^2}$ barで他が1 barでも平衡状態になる。他にも、いくらでも平衡状態を与えるガス分圧の組み合わせを考えることができて、その数は無数にある。つまり、平衡定数の値がひとつに決まっても、平衡ガス分圧は一意に決まるわけではない。

さて、我々は何を求めたかったか、取り扱いたかった問題に戻ってみよう。反応管に入れるメタンと水の流量から、反応が平衡まで達した時の各ガスの流量が知りたかったわけである。「反応管に入れる流量」が与えられているということは、単位時間当りに反応管に入ってくる炭素、水素、酸素原子の数が決まっているということであり、質量保存の法則から平衡に達した後もこの数は維持されなければならない。つまり、炭素、水素、酸素、それぞれについての3つの質量保存の式が、式（3）に加わる。これにさらに、全圧がいくらかという条件が加わることで、平衡状態のガス分圧が一意に決まる、というわけである。

ここまでをまとめてみよう。平衡状態のメタン転化率は、

平衡定数と平衡ガス分圧の釣り合い式（式（3））
ガス流入量と平衡ガス流出量との間の質量保存則
系の全圧

から、求めることになる。

メタン水蒸気改質反応以外の反応

メタン水蒸気改質反応に伴って、水性シフト反応や炭素析出反応が生じる。これらを考慮して、平衡状態を計算しなければならない。原理的には、考えうる全ての反応を考慮しなければならないが、それは大変なので、あまり起こらない反応（メタノールが生成する逆メタノール分解反応など）は考慮しない。また、熱力学的には十分起こりうるが、扱っている触媒上ではあまり起こらない反応についても、これを考慮しない場合の平衡と比べた方が、触媒の性能を評価しやすい場合がある。ここでは、我々が現在扱っている触媒実験の結果に即して、炭素析出が起こらない場合の平衡を考えることにする（ただし、必要ならば、炭素析出が起こる場合の計算は、CEAパッケージでできる）。

また、既に考慮に入れている反応から演繹できるようなものも、考える必要がない。例えば、メタンと二酸化炭素から水素を生成するメタンドライ改質反応は、メタン水蒸気改質反応と逆水性シフト反応が同時に起こったのと同じなので、この二つの反応を考慮しておけば、十分である。

従って、ここでは、メタン水蒸気改質反応以外の反応として、水性シフト反応だけを考える。反応式と平衡定数は、以下で与えられる。

$\begin{displaymath} \ensuremath{\xi_{\mathrm{W}}}: \ensuremath{\mathrm{CO}}+ \en... ...hrm{O}}= \ensuremath{\mathrm{CO}_2}+ \ensuremath{\mathrm{H}_2} \end{displaymath}$

(4)

$\begin{displaymath} \ensuremath{K_{\mathrm{W}}}= \frac{P_{\ensuremath{\mathrm{CO... ...h{\mathrm{CO}}} \cdot P_{\ensuremath{\mathrm{H}_2\mathrm{O}}}} \end{displaymath}$

(5)

反応進行度 $\xi$

反応進行度 $\xi$ という変数を導入する。既に、反応式（1）及び（4）の中に記号だけ書いておいた。 $\xi$ を求めるべき変数とする方法は、自然に質量保存則を満足させながら、扱う方程式の数を減らせて、便利である。

反応進行度 $\xi$ は、どれだけ反応が進行したかをモル単位で表したものである。例えば、式（1）のメタン水蒸気改質反応で考えてみよう。反応進行度 $\ensuremath{\xi_{\mathrm{M}}}$ が1 molということは、反応管の中から1 molのメタンと水が消えて、1 molの一酸化炭素と3 molの水素が生成するということを意味する（そのように決めるわけだ）。反応式は、左辺と右辺で各元素の数が変わらないようにできているので、反応進行度に応じて各ガスのモル数を変化させても、必ず、質量が保存される。反応進行度を導入することで、質量保存則を保った上で、一気に独立変数を各反応式につき、一つに減らせるわけである。

モル数で記述する

モル数で平衡状態を記述してみよう。各ガスのモル数

とガス分圧

は、以下の状態方程式に従う。

$\begin{displaymath} P_i V = n_iRT \; (i = \ensuremath{\mathrm{CH}_4}, \ensuremat... ...thrm{CO}_2}, \ensuremath{\mathrm{H}_2}, \mathrm{and so on}) \end{displaymath}$

(6)

この平衡状態でのモル数

を求めたいわけである。

式（6）を、式（3）と式（5）に代入すると、それぞれ、

$\begin{displaymath} \ensuremath{K_{\mathrm{M}}}= \frac{n_{\ensuremath{\mathrm{CO... ...th{\mathrm{H}_2\mathrm{O}}}} \cdot \left(\frac{RT}{V}\right)^2 \end{displaymath}$

(7)

$\begin{displaymath} \ensuremath{K_{\mathrm{W}}}= \frac{n_{\ensuremath{\mathrm{CO... ...h{\mathrm{CO}}} \cdot n_{\ensuremath{\mathrm{H}_2\mathrm{O}}}} \end{displaymath}$

(8)

となる。メタン水蒸気改質反応のように、反応後に系のモル数が増減する場合¹は、

という項がつく。これが、後々、定圧条件（流通系）と定体積条件（密封系）とで、計算方法を分ける原因となる。

さて、これを初期のモル数を所与として、反応進行度を変数とする式に直そう。平衡状態のモル数を初期のモル数と反応進行度 $\xi$ で表すと、次のようになる。

$\displaystyle \begin{array}{lll} \jot=20pt n_{\ensuremath{\mathrm{CH}_4}} &=& n... ... n^0_{\ensuremath{\mathrm{CO}_2}} + {\ensuremath{\xi_{\mathrm{W}}}} \end{array}$

(9)

反応が $\ensuremath{\xi_{\mathrm{M}}}$ 及び $\ensuremath{\xi_{\mathrm{W}}}$ だけ進むと、どのガスがどれだけ増えて、どれだけ減るか、二つの反応式（1）、（4）をじっと見ると、上の関係は、了解できるだろう。全体のモル量

は、初期の全モル量

とメタン水蒸気改質反応の反応進行度 $\ensuremath{\xi_{\mathrm{M}}}$ から、

$\begin{displaymath} n = n^0 + 2 \ensuremath{\xi_{\mathrm{M}}} \end{displaymath}$

(10)

となる。

式（9）を、平衡定数の式（7）、( 8) に代入すると、

$\begin{displaymath} \ensuremath{K_{\mathrm{M}}}= \frac{(n^0_{\ensuremath{\mathrm... ...uremath{\xi_{\mathrm{W}}}})} \cdot \left(\frac{RT}{V}\right)^2 \end{displaymath}$

(11)

$\begin{displaymath} \ensuremath{K_{\mathrm{W}}}= \frac{(n^0_{\ensuremath{\mathrm... ...uremath{\xi_{\mathrm{M}}}} - {\ensuremath{\xi_{\mathrm{W}}}})} \end{displaymath}$

(12)

となる。

この2式を満足するように、2つの反応進行度を解けば良い。ごちゃごちゃしてきたが、反応進行度との他は、初期値だったり、所与の条件、定数として、与えられている。

さらにモル比へ

流通系の場合、流速を2倍にしても、系の圧力が変わらない限り、平衡状態（出口）でのガス組成は変わらない。密封系で定容系反応の組み合わせでも、モル量を2倍にしても、圧力の影響を受けず、平衡組成は変わらない。密封系で、かつ、非定容系反応では、投入するガスのモル量を2倍にすると、圧力があがる影響を受けて、平衡状態でのガス組成は変わる。例えば、メタン水蒸気改質反応は起こりにくくなる。しかし、これは初期の全圧を高くした影響として表現できる。何が言いたいかというと、流通系でも密封系でも、定容系でも非定容系でも、（つまりどんな場合でも）初期のガス組成と全圧を与えれば、平衡状態のガス組成と全圧は決まるということである。モル量で初期値を与える必要はないのだ。ガス組成、つまり、モル比で扱う方が、モル数で扱うよりも、実験結果との比較が簡単にできる。メタンと水蒸気の比が分かっていれば、25 ml/minが何mol/minなのかをわざわざ計算しなくても良いというのは、なかなか便利なのだ。

具体的には、これまでモル量で表現してきた式（11）及び式（12）を、初期のモル比（は初期全モル量）で表現する。すると、以下のようになる。

$\begin{displaymath} \ensuremath{K_{\mathrm{M}}}= \frac{(x^0_{\ensuremath{\mathrm... ...emath{y_{\mathrm{W}}}})} \cdot \left(\frac{n^0 RT}{V}\right)^2 \end{displaymath}$

(13)

$\begin{displaymath} \ensuremath{K_{\mathrm{W}}}= \frac{(x^0_{\ensuremath{\mathrm... ...\ensuremath{y_{\mathrm{M}}}} - {\ensuremath{y_{\mathrm{W}}}})} \end{displaymath}$

(14)

ここで、 $\ensuremath{y_{\mathrm{M}}}$ と $\ensuremath{y_{\mathrm{W}}}$ は、反応の進み具合を、初期の全モル量で規格化したもので、

$\displaystyle \begin{array}{lll} \ensuremath{y_{\mathrm{M}}}&=& \ensuremath{\xi... ... \ensuremath{y_{\mathrm{W}}}&=& \ensuremath{\xi_{\mathrm{W}}}/ n^0 \end{array}$

(15)

である。この二つの無次元数を求めることになる。

定容系反応である水性シフト反応は、分母と分子のモル数にかかかる次数の合計が一致するので、が式から消える。非定容系反応であるメタン水蒸気改質反応は、そういうわけにはいかず、が残る。この項目の扱いは、定体積条件の方が簡単なので、そちらから説明を始める。

定体積条件

式（13）の右辺に残った

は、定積条件では

（

は初期体積）と書くことができて、所与の定数となる。理想気体の状態方程式から、

$\begin{displaymath} \left( \frac{n^0 RT}{V^0} \right)^2 = {P^0}^2 \end{displaymath}$

(16)

と、一気に初期の圧力

で置き換えよう。これを式（13）に代入して、

$\begin{displaymath} \ensuremath{K_{\mathrm{M}}}= \frac{(x^0_{\ensuremath{\mathrm... ...h{y_{\mathrm{M}}}- \ensuremath{y_{\mathrm{W}}})} \cdot {P^0}^2 \end{displaymath}$

(17)

を得る。

定体積条件は、密封系の実験に対応する。初期圧力が高いと、式（17）から、生成物（式の分子）の積が小さくなることが分かる。つまりは、反応が進行しにくいということであり、直感ともあう。

定圧条件

定圧条件では、式（13）右辺の

は、変数である。この中の

がもし

だったら、以下のように（平衡状態が達成された状態での）気体の状態方程式が適用できる。

$\begin{displaymath} \left( \frac{n RT}{V} \right)^2 = {P^0}^2 \end{displaymath}$

(18)

ここで、右辺を

としたのは、定圧条件

だからである。

式（18）をの置き換えに利用するために、を $n=n^0 \cdot n/n^0$ と考えてみよう。さらに、と定義して（全モル量変化と呼び）、 $n=n^0 \cdot m$ とする。すると、式（18）は、

$\begin{displaymath} \left( \frac{n^0 m RT}{V} \right)^2 = {P^0}^2 \end{displaymath}$

(19)

となり、これから

を次のように置き換えることができる。

$\begin{displaymath} \left( \frac{n^0 RT}{V} \right)^2 = \left(\frac{P^0}{m}\right)^2 \end{displaymath}$

(20)

これを式（13）に代入して、

$\begin{displaymath} \ensuremath{K_{\mathrm{M}}}= \frac{(x^0_{\ensuremath{\mathrm... ...uremath{y_{\mathrm{W}}})} \cdot \left( \frac{P^0}{m} \right)^2 \end{displaymath}$

(21)

を得る。

最後に、全モル量変化を書き下しておこう。 $n=n^0 + 2 \ensuremath{\xi_{\mathrm{M}}}$ （式（10））の両辺をで割って、

$\begin{displaymath} m = \frac{n}{n^0} = 1 + 2 \ensuremath{y_{\mathrm{M}}} \end{displaymath}$

(22)

となる。従って、最終的に以下の式を得る。

$\begin{displaymath} \ensuremath{K_{\mathrm{M}}}= \frac{(x^0_{\ensuremath{\mathrm... ...\left( \frac{P^0}{1 + 2 \ensuremath{y_{\mathrm{M}}}} \right)^2 \end{displaymath}$

(23)

定体積条件と比べると、 $1 + 2 \ensuremath{y_{\mathrm{M}}}$ の分だけ、 $\ensuremath{\mathrm{CO}}$ や $\ensuremath{\mathrm{H}_2}$ が多く生成することになる。全モル量が増えるメタン水蒸気改質反応は、定圧条件、つまり、体積膨張が許されている場合の方が、より反応が進行するということであり、これは直感ともあう。定圧条件に対応する実験は、例えば、流通系である。本稿で計算したい流通系は、従って、この式を用いることになる。

解き方

メタン水蒸気改質反応の平衡状態を表す式（17）（定積条件）あるいは（23）（定圧条件）と、水性シフト反応を表す式（14）を連立させて、反応進行度を全モル量で規格化した $\ensuremath{y_{\mathrm{M}}}$ と $\ensuremath{y_{\mathrm{W}}}$ を求めれば良い。いろいろやってみたが、4次方程式と2次方程式の連立となり、陽に解を提示することはできなかった。ニュートンラプソン法など、数値解析で近似解を求めることになる。ここでその詳細は示さないが、解き方の方針と注意点を述べておく。

まず、式（14）を $\ensuremath{y_{\mathrm{W}}}$ の2次方程式として解く。これを式（17）あるいは（23）に代入して、 $\ensuremath{y_{\mathrm{M}}}$ だけを変数とする方程式をたてる。これに対して、ニュートンラプソン法を適用するわけである。

注意点は、（1）解として不適格なものを排除することと、（2）ニュートンラプソン法を適用する式の形態や初期値を適切に設定することである。（1）を考えてみよう。解は複数存在するが、そのうち現実に即しているものは一つだけのはずである。現実に即さない解を排除しなければならない。式（9）で、左辺の各物質のモル量が、負の値になることはない。これを両辺で割ると、 $\ensuremath{y_{\mathrm{M}}}$ や $\ensuremath{y_{\mathrm{W}}}$ の解の範囲が分かる。例えば、これらはそれぞれ、初期のメタンや水蒸気のモル比を超えてはいけない。他に初期条件から常識的に判断できることもある。初期にメタンと水蒸気を流す場合、メタン水蒸気改質反応や水性シフト反応が「進行」すると期待でき、決して、「後退」＝「逆反応」は起こらないだろう。つまりは、 $\ensuremath{y_{\mathrm{M}}}$ や $\ensuremath{y_{\mathrm{W}}}$ は正の値を取るはずである。こういうところから、解を限定していくわけである。

式（14）を $\ensuremath{y_{\mathrm{W}}}$ の2次方程式として解いたときの、適切な解の選び方を紹介しておこう。2次方程式の解の方程式から、以下の2組の解を得る。

$\begin{displaymath} \ensuremath{y_{\mathrm{W}}} =\frac{A \pm \sqrt{B}} {2({\it\ensuremath{K_{\mathrm{W}}}}-1)} \end{displaymath}$

(24)

$\displaystyle \begin{array}{l} A = \ensuremath{K_{\mathrm{W}}}x^0_{\ensuremath{... ...nsuremath{K_{\mathrm{W}}}+9 \right) {\ensuremath{y_{\mathrm{M}}}}^2 \end{array}$

(25)

式（24）のどちらの解を選ぶか。本稿で考えているような、初期にメタンと水蒸気を反応物とするケースでは、 $\ensuremath{y_{\mathrm{W}}}$ が負となることはない。 $\ensuremath{K_{\mathrm{W}}}$ の値を見ると、600では $\ensuremath{K_{\mathrm{W}}}\simeq 0.97$ で1より小さく、温度の上昇とともに値は小さくなる。式（24）の分母は従って、600以上では負になる。 $\ensuremath{y_{\mathrm{W}}}$ を正にするには、式（24）の分子は負でなければならないだろう。も $\sqrt{B}$ も複雑ではあるがどちらも正の数なので、足すと正になり、これはあり得ない。とすると、取りうる解は、以下ということになる。

$\begin{displaymath} \ensuremath{y_{\mathrm{W}}} =\frac{A - \sqrt{B}} {2({\it\ensuremath{K_{\mathrm{W}}}}-1)} \end{displaymath}$

(26)

これを式（17）あるいは（23）に代入して、ニュートンラプソン法を適用することになるわけである。この際の注意点は、式の形である。式（17）あるいは（23）の分母には変数を含む項がある。 $\ensuremath{K_{\mathrm{M}}}$ は1より大きいので、この項がかなり小さくなってしまい、上手く集束しない。分母の、変数を含む項をあらかじめ両辺にかけて分母から消しておけば、上手く集束して、解が求まる。初期値は、あり得る範囲で取れば良い。筆者は、初期値にを用いており、それで上手く計算できている。

実勢の計算は、方程式を解くプログラムを利用すると良い。筆者は、「maxima」というオープンソースのプログラムを利用した²。その際のコードを付録（7.1節）につけておく。

平衡定数の温度依存性

平衡定数は、エンタルピーや比熱などの熱力学データから計算できる。温度の何次までの関数とするかで、式や温度定数の値が変わる。ここでは、有限会社コムテック・クウェストが提供している資料³を参考にした。

自由エネルギーを最小化する方法

自由エネルギーを最小化する方法については、松本・横川らの日本語の解説に詳しい⁴。また、CEAパッケージの著者らによる解説も、参考になるだろう⁵。

以下、自由エネルギーを最小化する方法で計算するCEAパッケージについて解説する。

CEAパッケージ

CEAパッケージは、NASAのBonnie McBride博士とSanford Gordon博士のグループが開発した平衡計算用のコードである⁶。豊富な熱力学データが用意されていて、かなりのガス反応を扱うことができる。コードは、JAVAで書かれており、Windows、Mac、Unixで使用できる。Windows版のインストールと使い方は、東京大学准教授三好明先生のwebページ⁷に詳しい。Mac版も同じようにインストールできるが、GUI版の実行ファイルは上手く動かないので、ソースからコンパイルする必要がある。

使い方は、先に紹介した三好先生の解説が分かりやすいだろう。

メタン水蒸気改質反応について計算した時の、入力ファイルを参考までに載せておく。

1. problem   case=SC1
2.   tp   t,c=500,550,600,650,700,750,800,850,900,  p,bar=1,
3. react  
4.   name=CH4 moles=1  
5.   name=H2O moles=1  
6. only 
7.   CH4 CO CO2 H2 H2O 
8. output  short 
9. end

1行目は、この計算の名前で、ここは何でも良い。2行目の"tp"は、温度と圧力を一定にするという計算条件を指示している。その後に、計算する温度と圧力を設定している。もしも定体積条件で計算したい場合は、"tv"と書き、その後、温度と体積を指定するわけである。

3行目から反応物を指定している。ここでは、ガスの名前とモル数で指定している。これは、の例である。

6行目が大事である。ここでは、生成物として考慮するガスの種類を指定してある。これを指定しなくても計算でき、その方がよりエネルギーの低い熱力学的安定状態を与える。しかし、2.3説で述べたように、起こる反応を限定して、平衡を計算したい場合がある。その時に、この"only"を利用しよう。ここでは、メタン水蒸気改質反応と水性シフト反応で現れるガスだけに計算対象を限定している。

8行目は出力の形式の指定で、9行目は終了コードである。

このように、反応条件を指定して、生成物を与え、必要に応じて反応を限定するだけで、平衡が計算できる有り難いプログラムである。

比較

**図1:** 平衡定数から計算した平衡メタン転化率（実線、破線、一点鎖線）とCEAパッケージを利用した計算値（丸、四角、三角）の比較

平衡定数から計算した結果とCEAで計算した結果を比較してみよう。図1は、種々の比、全圧における平衡メタン転化率を計算したものである。炭素の析出はないと仮定している。また、窒素ガスなどのバッファーガスは入れていない。平衡定数からの計算結果（実線、破線、一点鎖線）は、CEAを利用した自由エネルギーを最小化する計算結果（丸、四角、三角）とよく一致している。ここで計算した値は、Rostrup-Nielsen⁸が報じている値とも良く一致した。ただし、彼は計算の方法や前提（炭素析出やバッファーガスの有無）を記していないので、ここで詳しく比較の結果を述べるのはやめる。

参考にした文献

参考文献は、適時、脚注に載せた。それ以外の、考え方を教えてくれた、重要な2つの文献を挙げておく。

Cpt. 6, "Physical Chemistry" by Walter J. Moore.
「化学ポテンシャルと平衡定数」小口達夫、梶本興亜、山崎勝義著⁹

付録

maximaコード

メタンと水の比、初期の全圧を与えて、指定した温度範囲でメタン転化率を計算するmaximaコードを載せる。定圧条件で計算している。定体積条件の式は、コメントアウトしてある。

/*newton*/
load(newton);

/*MSR equilibrium constant:*/
define(lnKm(T),-delH0m/R/T + delAm/R*log(T) + delBm/2/R*T + delCm/6/R*T^2 + 
delDm/12/R * T^3 + delEm/20/R*T^4+Im);
delH0m:1.93E5;
delAm:36.878;
delBm:1.02E-1;
delCm:-3.17E-4;
delDm:2.54E-7;
delEm:-6.07E-11;
R:8.314;
Im:-6.1697;

/*WGS equilibrium constant:*/
define(lnKw(T),-delH0w/R/T + delAw/R*log(T) + delBw/2/R*T + delCw/6/R*T^2 + 
delDw/12/R * T^3 + delEw/20/R*T^4+Iw);}
delH0w:-4.06E4;
delAw:-10.653;
delBw:7.75E-2;
delCw:-1.08E-4;
delDw:6.59E-8;
delEw:-1.50E-11;
Iw:1.2438;

/*definition used*/
CH4:CH4_0-Ym;
H2O:H2O_0-Ym-Yw;
H2:H2_0+3*Ym+Yw;
CO:CO_0+Ym-Yw;
CO2:CO2_0+Yw;

/*MSR equilibrium equation:*/
MSR:Km*CH4*H2O*(1+2*Ym)^2 - H2^3*CO*P0^2;/*Constant pressure*/
/*MSR:Km*CH4*H2O - H2^3*CO*P0^2;Constant volume*/

/*WGS equilibrium equation:*/
WGS:Kw*CO*H2O - H2*CO2;

/*WGS solution*/
solWGS:solve(WGS,Yw);

/*Insert the solution of Yw into MSR, answer #1 should give a reasonable solution*/
rMSR:ev(MSR,Yw:rhs(solWGS[1]));

/*Calculate conversion as a function of Tempratre*/
MSRvsT(Ts,Te,Tstep,SoverC,pressure):=block([conversion,Tstart:Ts+273.15,
Tend:Te+273.15,tempMSR,guess:0,steam,methane,tempANS],
/*methane ans steam*/
tempANS:solve([s+c=1,s/c=SoverC],[s,c])[1],
steam:rhs(tempANS[1]),
methane:rhs(tempANS[2]),
for i:Tstart thru Tend step Tstep do(
tempMSR:ev(rMSR,CH4_0:methane,H2O_0:steam,P0:pressure,CO_0:0,H2_0:0,CO2_0:0,
Km:exp(lnKm(i)),Kw:exp(lnKw(i)),T:i),
solYm:newton(tempMSR,guess),
/*solYw:rhs(ev(solWGS[1],Ym:solYm,H2O_0:steam,CO_0:0,H2_0:0,CO2_0:0,
Km:exp(lnKm(i)),Kw:exp(lnKw(i)),T:i)),
print("Ymsr = ",solYm, ", Ywgs = ", solYw),*/
/*Methane conversion*/
conversion:solYm/methane*100,
print(round(i-273.15),float(round(conversion*10)/10))),
return("done."));

MSRvsT(600,900,10,2,1);

この文書について...

メタン水蒸気改質反応の平衡計算

この文書はLaTeX2HTML 翻訳プログラム Version 2002-2-1 (1.71)

を日本語化したもの( 2002-2-1 (1.71) JA patch-1.9 版)

を用いて生成されました。

コマンド行は以下の通りでした。:
latex2html -no_navigation -split +0 MSR-equilibrium-demura.

翻訳は Masahiko Demura によって平成24年1月13日に実行されました。

脚注

... となる。メタン水蒸気改質反応のように、反応後に系のモル数が増減する場合 ¹: 反応で系のモル数が変化しない場合を定容系、変化する場合を非定容系という。
...実勢の計算は、方程式を解くプログラムを利用すると良い。筆者は、「maxima」というオープンソースのプログラムを利用した ²: URL http://maxima.sourceforge.net/
... 平衡定数は、エンタルピーや比熱などの熱力学データから計算できる。温度の何次までの関数とするかで、式や温度定数の値が変わる。ここでは、有限会社コムテック・クウェストが提供している資料 ³: URL http://comtecquest.com/_userdata/EqulibriumConstant&SteamReforming.pdf
...自由エネルギーを最小化する方法については、松本・横川らの日本語の解説に詳しい ⁴: T. Matsumoto, H. Yokokawa, Netsu Sokutei 19 (1992) 170173.
...。また、CEAパッケージの著者らによる解説も、参考になるだろう ⁵: S. Gordon, B.J. McBridge, Computer Program for Calculation of Complex Equilibrium Compositions and Applications: Part I. Analysis, NASA Reference Publication, 1994.
... Gordon博士のグループが開発した平衡計算用のコードである ⁶: URL http://www.grc.nasa.gov/WWW/CEAWeb/
...。豊富な熱力学データが用意されていて、かなりのガス反応を扱うことができる。コードは、JAVAで書かれており、Windows、Mac、Unixで使用できる。Windows版のインストールと使い方は、東京大学准教授三好明先生のwebページ ⁷: http://www.frad.t.u-tokyo.ac.jp/public/cea/cea_index.html.ja
...における平衡メタン転化率を計算したものである。炭素の析出はないと仮定している。また、窒素ガスなどのバッファーガスは入れていない。平衡定数からの計算結果（実線、破線、一点鎖線）は、CEAを利用した自由エネルギーを最小化する計算結果（丸、四角、三角）とよく一致している。ここで計算した値は、Rostrup-Nielsen ⁸: J. Rostrup-Nielsen, in:, Phys Chem Chem Phys, 2001, pp. 283288.
... 著 ⁹: ダウンロードURL http://home.hiroshima-u.ac.jp/kyam/pages/results/monograph/Refthermo_G.pdf

Masahiko Demura 平成24年1月13日