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

メタン水蒸気改質反応の平衡計算
出村雅彦
2011 年 12 月 22 日
目次
1
取り扱う問題とアプローチ
1
2
平衡定数から計算する方法
2
2.1
平衡定数と平衡ガス分圧 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2.2
平衡状態を求めるということ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2.3
メタン水蒸気改質反応以外の反応 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.4
反応進行度 ξ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.5
モル数で記述する . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.6
さらにモル比へ
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.7
定体積条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.8
定圧条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.9
解き方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
3
平衡定数の温度依存性
8
4
自由エネルギーを最小化する方法
8
CEA パッケージ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
5
比較
9
6
参考にした文献
10
7
付録
10
maxima コード . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
4.1
7.1
1 取り扱う問題とアプローチ
反応管にメタンを 25 ml/min、水蒸気を 25 ml/min 流そう。反応管を加熱するとメタン水蒸気改質反応や
水性シフト反応が起こる。化学的平衡に達すると、出口での各ガスの流量、そして、メタン転化率はどうなる
か？—これが、ここで取り扱う問題である。メタン・水蒸気比（S/C）や不活性バッファーガス（N2 等）導入
1
が、平衡メタン転化率にどんな影響を及ぼすのかも、視野のうちに入れたい。
平衡状態では、それ以上反応が進行せず、Gibbs エネルギーが極小化する。Gibbs エネルギーが極小化する
反応の進行度合いを求めれば良いわけだ。その方法には、二つある。一つは、平衡定数から求める方法。これ
は、昔からある古典的な方法である。ここでは、この方法について、詳しく述べる。もう一つの方法は、Gibbs
エネルギーを直接最小化する方法である。これは平衡の定義そのものである。この方法は、フリーで公開され
ているパッケージ（CEA）を利用すると簡単に試すことができる。ここでは、その使い方を紹介するだけに
とどめ、詳しい解説はしない（というか、勉強不足でできない）
。もちろん、基本的にどちらの方法も、同じ答
えを与える。
では、まずは、平衡定数から計算する方法から、見てみよう。古典な方法ではあるが、これを理解しておく
と、計算結果をより深く理解できると思う。
2 平衡定数から計算する方法
2.1 平衡定数と平衡ガス分圧
ある反応の平衡定数は、反応式右辺（反応物）のガスの平衡分圧を分母に、左辺（生成物）の平衡ガス分圧
を分母においた形をしている。メタン水蒸気改質反応で説明しよう。反応式は、以下で表される。
ξM : CH4 + H2 O = CO + 3H2
(1)
ξM は反応の進行度を表す変数で、後ほど登場する。この反応式で書いた時のメタン水蒸気改質反応の平衡定
数は、
KM =
KM =
PCO
P0
·
PCH4
P0
(
PH2
P0
)3
PH2 O
P0
PH2 3
·
PCO ·
PCH4 · PH2 O
(2)
(3)
となる。ここで、Pi （i は物質名）は、平衡に達した時の部分モル分圧（以下、ガス分圧と呼ぶ）、P 0 は標準
状態の圧力で通常は 1 bar である。二つ式があるが、P 0 = 1 bar なので、計算する上で違いはない。以下で
は、式（3）を採用する。ただし、各ガスのガス分圧は 1 bar に対する比で表されていると考えて、KM は無
次元数であることを忘れないようにしたい。
水素のガス分圧が 3 乗となっていることに注意してほしい。反応式中の係数でベキ乗するわけである。とい
うことは、式（1）を 2 倍すると、その平衡定数は式（3）を 2 乗したものになる。つまり、反応式を与えない
と平衡定数は決まらない。メタン水蒸気改質反応の場合はわざわざ 2 倍した式を使うことはないが、一酸化炭
素の酸化反応（CO シフト反応）などでは、一酸化炭素の係数を 1 にして書いてみたり、分数が出てくるのを
避けるために係数を 2 にしてみたり、目的によって反応式も変わり、それに伴って、平衡定数の値も変わるこ
とになるので、注意が必要である。
2.2 平衡状態を求めるということ
平衡定数は、
（純粋な）各ガス分子の標準状態（1 bar）の Gibbs 自由エネルギー（化学ポテンシャル）から
求めることができ、反応と温度を決めれば、一意にその値が決まる。その値に合うように、各ガス分子のガス
2
分圧を決めれば良いわけである。これが決まれば、平衡状態の流量が計算でき、メタンの平衡転化率や各ガス
の平衡選択率が計算できる。
例えば、800˚C のメタン水蒸気改質反応の平衡定数は、KM ≃ 1.57 · 102 である。この値になるように、ガス
分圧を決めてみよう。PCO = 1.57 · 102 bar で、他が全て 1 bar というガス分圧の組み合わせは、平衡定数と
一致する。従って、ひとつの平衡状態である。似たような解として、PH2 =
√
3
1.57 · 102 bar で他が 1 bar で
も平衡状態になる。他にも、いくらでも平衡状態を与えるガス分圧の組み合わせを考えることができて、その
数は無数にある。つまり、平衡定数の値がひとつに決まっても、平衡ガス分圧は一意に決まるわけではない。
さて、我々は何を求めたかったか、取り扱いたかった問題に戻ってみよう。反応管に入れるメタンと水の流
量から、反応が平衡まで達した時の各ガスの流量が知りたかったわけである。「反応管に入れる流量」が与え
られているということは、単位時間当りに反応管に入ってくる炭素、水素、酸素原子の数が決まっているとい
うことであり、質量保存の法則から平衡に達した後もこの数は維持されなければならない。つまり、炭素、水
素、酸素、それぞれについての 3 つの質量保存の式が、式（3）に加わる。これにさらに、全圧がいくらかと
いう条件が加わることで、平衡状態のガス分圧が一意に決まる、というわけである。
ここまでをまとめてみよう。平衡状態のメタン転化率は、
1. 平衡定数と平衡ガス分圧の釣り合い式（式（3））
2. ガス流入量と平衡ガス流出量との間の質量保存則
3. 系の全圧
から、求めることになる。
2.3 メタン水蒸気改質反応以外の反応
メタン水蒸気改質反応に伴って、水性シフト反応や炭素析出反応が生じる。これらを考慮して、平衡状態を
計算しなければならない。原理的には、考えうる全ての反応を考慮しなければならないが、それは大変なの
で、あまり起こらない反応（メタノールが生成する逆メタノール分解反応など）は考慮しない。また、熱力学
的には十分起こりうるが、扱っている触媒上ではあまり起こらない反応についても、これを考慮しない場合の
平衡と比べた方が、触媒の性能を評価しやすい場合がある。ここでは、我々が現在扱っている触媒実験の結果
に即して、炭素析出が起こらない場合の平衡を考えることにする（ただし、必要ならば、炭素析出が起こる場
合の計算は、CEA パッケージでできる）。
また、既に考慮に入れている反応から演繹できるようなものも、考える必要がない。例えば、メタンと二酸
化炭素から水素を生成するメタンドライ改質反応は、メタン水蒸気改質反応と逆水性シフト反応が同時に起
こったのと同じなので、この二つの反応を考慮しておけば、十分である。
従って、ここでは、メタン水蒸気改質反応以外の反応として、水性シフト反応だけを考える。反応式と平衡
定数は、以下で与えられる。
ξW : CO + H2 O = CO2 + H2
KW =
PCO2 · PH2
PCO · PH2 O
3
(4)
(5)
2.4 反応進行度 ξ
反応進行度 ξ という変数を導入する。既に、反応式（1）及び（4）の中に記号だけ書いておいた。ξ を求め
るべき変数とする方法は、自然に質量保存則を満足させながら、扱う方程式の数を減らせて、便利である。
反応進行度 ξ は、どれだけ反応が進行したかをモル単位で表したものである。例えば、式（1）のメタン水
蒸気改質反応で考えてみよう。反応進行度 ξM が 1 mol ということは、反応管の中から 1 mol のメタンと水
が消えて、1 mol の一酸化炭素と 3 mol の水素が生成するということを意味する（そのように決めるわけだ）。
反応式は、左辺と右辺で各元素の数が変わらないようにできているので、反応進行度に応じて各ガスのモル数
を変化させても、必ず、質量が保存される。反応進行度を導入することで、質量保存則を保った上で、一気に
独立変数を各反応式につき、一つに減らせるわけである。
2.5 モル数で記述する
モル数で平衡状態を記述してみよう。各ガスのモル数 ni とガス分圧 Pi は、以下の状態方程式に従う。
Pi V = ni RT (i = CH4 , H2 O, CO, CO2 , H2 , and so on)
(6)
この平衡状態でのモル数 ni を求めたいわけである。
式（6）を、式（3）と式（5）に代入すると、それぞれ、
KM
nCO · nH2 3
=
·
nCH4 · nH2 O
KW =
(
RT
V
)2
nCO2 · nH2
nCO · nH2 O
(7)
(8)
となる。メタン水蒸気改質反応のように、反応後に系のモル数が増減する場合*1 は、RT /V という項がつく。
これが、後々、定圧条件（流通系）と定体積条件（密封系）とで、計算方法を分ける原因となる。
さて、これを初期のモル数を所与として、反応進行度を変数とする式に直そう。平衡状態のモル数 ni を初
期のモル数 n0i と反応進行度 ξ で表すと、次のようになる。
nCH4
n H2 O
n H2
nCO
nCO2
=
=
=
=
=
n0CH4 − ξM
n0H2 O − ξM − ξW
n0H2 + 3ξM + ξW
n0CO + ξM − ξW
n0CO2 + ξW
(9)
反応が ξM 及び ξW だけ進むと、どのガスがどれだけ増えて、どれだけ減るか、二つの反応式（1）
、
（4）をじっ
と見ると、上の関係は、了解できるだろう。全体のモル量 n は、初期の全モル量 n0 とメタン水蒸気改質反応
の反応進行度 ξM から、
n = n0 + 2ξM
となる。
*1
反応で系のモル数が変化しない場合を定容系、変化する場合を非定容系という。
4
(10)
式（9）を、平衡定数の式（7）、( 8) に代入すると、
3
KM
(n0 + ξM − ξW ) · (n0H2 + 3ξM + ξW )
= CO 0
·
(nCH4 − ξM ) · (n0H2 O − ξM − ξW )
KW =
(
RT
V
)2
(11)
(n0CO2 + ξW ) · (n0H2 + 3ξM + ξW )
0
(nCO + ξM − ξW ) · (n0H2 O − ξM − ξW )
(12)
となる。
この 2 式を満足するように、2 つの反応進行度を解けば良い。ごちゃごちゃしてきたが、反応進行度と V の
他は、初期値だったり、所与の条件、定数として、与えられている。
2.6 さらにモル比へ
流通系の場合、流速を 2 倍にしても、系の圧力が変わらない限り、平衡状態（出口）でのガス組成は変わら
ない。密封系で定容系反応の組み合わせでも、モル量を 2 倍にしても、圧力の影響を受けず、平衡組成は変わ
らない。密封系で、かつ、非定容系反応では、投入するガスのモル量を 2 倍にすると、圧力があがる影響を受
けて、平衡状態でのガス組成は変わる。例えば、メタン水蒸気改質反応は起こりにくくなる。しかし、これは
初期の全圧を高くした影響として表現できる。何が言いたいかというと、流通系でも密封系でも、定容系でも
非定容系でも、（つまりどんな場合でも）初期のガス組成と全圧を与えれば、平衡状態のガス組成と全圧は決
まるということである。モル量で初期値を与える必要はないのだ。ガス組成、つまり、モル比で扱う方が、モ
ル数で扱うよりも、実験結果との比較が簡単にできる。メタンと水蒸気の比が分かっていれば、25 ml/min が
何 mol/min なのかをわざわざ計算しなくても良いというのは、なかなか便利なのだ。
具体的には、これまでモル量 ni で表現してきた式（11）及び式（12）を、初期のモル比 x0i = n0i /n（n0 は
初期全モル量）で表現する。すると、以下のようになる。
3
KM =
(x0CO + yM − yW ) · (x0H2 + 3yM + yW )
·
(x0CH4 − yM ) · (x0H2 O − yM − yW )
KW =
(
n0 RT
V
(x0CO2 + yW ) · (x0H2 + 3yM + yW )
(x0CO + yM − yW ) · (x0H2 O − yM − yW )
)2
(13)
(14)
ここで、yM と yW は、反応の進み具合を、初期の全モル量で規格化したもので、
yM
yW
= ξM /n0
= ξW /n0
(15)
である。この二つの無次元数を求めることになる。
定容系反応である水性シフト反応は、分母と分子のモル数にかかかる次数の合計が一致するので、n0 が式
から消える。非定容系反応であるメタン水蒸気改質反応は、そういうわけにはいかず、(n0 RT /V )2 が残る。
この項目の扱いは、定体積条件の方が簡単なので、そちらから説明を始める。
2.7 定体積条件
式（13）の右辺に残った (n0 RT /V )2 は、定積条件では (n0 RT /V 0 )2 （V 0 は初期体積）と書くことができ
て、所与の定数となる。理想気体の状態方程式から、
(
n0 RT
V0
)2
= P0
5
2
(16)
と、一気に初期の圧力 P 0 で置き換えよう。これを式（13）に代入して、
3
KM
(x0 + yM − yW ) · (x0H2 + 3yM + yW )
2
= CO 0
· P0
(xCH4 − yM ) · (x0H2 O − yM − yW )
(17)
を得る。
定体積条件は、密封系の実験に対応する。初期圧力 P 0 が高いと、式（17）から、 生成物（式の分子）の積
が小さくなることが分かる。つまりは、反応が進行しにくいということであり、直感ともあう。
2.8 定圧条件
定圧条件では、式（13）右辺の (n0 RT /V )2 は、変数である。この中の n0 がもし n だったら、以下のよう
に（平衡状態が達成された状態での）気体の状態方程式が適用できる。
(
nRT
V
)2
= P0
2
(18)
ここで、右辺を P 0 としたのは、定圧条件 P = P 0 だからである。
式（18）を (n0 RT /V )2 の置き換えに利用するために、n を n = n0 · n/n0 と考えてみよう。さらに、
m = n/n0 と定義して（全モル量変化と呼び）、n = n0 · m とする。すると、式（18）は、
(
n0 mRT
V
)2
= P0
2
(19)
となり、これから (n0 RT /V )2 を次のように置き換えることができる。
(
n0 RT
V
)2
(
=
P0
m
)2
(20)
これを式（13）に代入して、
3
KM =
(x0CO + yM − yW ) · (x0H2 + 3yM + yW )
·
(x0CH4 − yM ) · (x0H2 O − yM − yW )
(
P0
m
)2
(21)
を得る。
最後に、全モル量変化 m を書き下しておこう。n = n0 + 2ξM （式（10））の両辺を n0 で割って、
m=
n
= 1 + 2yM
n0
(22)
となる。従って、最終的に以下の式を得る。
3
KM =
(x0CO + yM − yW ) · (x0H2 + 3yM + yW )
·
(x0CH4 − yM ) · (x0H2 O − yM − yW )
(
P0
1 + 2yM
)2
(23)
定体積条件と比べると、1 + 2yM の分だけ、CO や H2 が多く生成することになる。全モル量が増えるメタ
ン水蒸気改質反応は、定圧条件、つまり、体積膨張が許されている場合の方が、より反応が進行するというこ
とであり、これは直感ともあう。定圧条件に対応する実験は、例えば、流通系である。本稿で計算したい流通
系は、従って、この式を用いることになる。
6
2.9 解き方
メタン水蒸気改質反応の平衡状態を表す式（17）
（定積条件）あるいは（23）
（定圧条件）と、水性シフト反
応を表す式（14）を連立させて、反応進行度を全モル量で規格化した yM と yW を求めれば良い。いろいろ
やってみたが、4 次方程式と 2 次方程式の連立となり、陽に解を提示することはできなかった。ニュートンラ
プソン法など、数値解析で近似解を求めることになる。ここでその詳細は示さないが、解き方の方針と注意点
を述べておく。
まず、式（14）を yW の 2 次方程式として解く。これを式（17）あるいは（23）に代入して、yM だけを変
数とする方程式をたてる。これに対して、ニュートンラプソン法を適用するわけである。
注意点は、
（1）解として不適格なものを排除することと、
（2）ニュートンラプソン法を適用する式の形態や
初期値を適切に設定することである。（1）を考えてみよう。解は複数存在するが、そのうち現実に即している
ものは一つだけのはずである。現実に即さない解を排除しなければならない。式（9）で、左辺の各物質のモ
ル量が、負の値になることはない。これを両辺 n0 で割ると、yM や yW の解の範囲が分かる。例えば、これら
はそれぞれ、初期のメタンや水蒸気のモル比を超えてはいけない。他に初期条件から常識的に判断できること
もある。初期にメタンと水蒸気を流す場合、メタン水蒸気改質反応や水性シフト反応が「進行」すると期待で
き、決して、
「後退」＝「逆反応」は起こらないだろう。つまりは、yM や yW は正の値を取るはずである。こ
ういうところから、解を限定していくわけである。
式（14）を yW の 2 次方程式として解いたときの、適切な解の選び方を紹介しておこう。2 次方程式の解の
方程式から、以下の 2 組の解を得る。
yW
√
A± B
=
2(KW − 1)
(24)
A = KW x0H2 O + x0H2 + KW x0CO + x0CO2 + 3yM
(
(
)
(
) )
2
B = KW 2 x0H2 O + 2KW x0H2 + 4KW − 2KW 2 x0CO + 2KW x0CO2 + 10KW − 4KW 2 yM x0H2 O
(
)
(
(
(25)
2
2
+x0H2 + 2KW x0CO + (4KW − 2) x0CO2 + 6yM x0H2 + KW 2 x0CO + 2KW x0CO2 + 4KW 2
(
)
2
+2KW ) yM ) x0CO + x0CO2 + (12KW − 6) yM x0CO2 + 4KW 2 − 4KW + 9 yM 2
式（24）のどちらの解を選ぶか。本稿で考えているような、初期にメタンと水蒸気を反応物とするケースで
は、yW が負となることはない。KW の値を見ると、600˚C では KW ≃ 0.97 で 1 より小さく、温度の上昇と
ともに値は小さくなる。式（24）の分母は従って、600˚C 以上では負になる。yW を正にするには、式（24）
の分子は負でなければならないだろう。A も
√
B も複雑ではあるがどちらも正の数なので、足すと正になり、
これはあり得ない。とすると、取りうる解は、以下ということになる。
yW =
√
A− B
2(KW − 1)
(26)
これを式（17）あるいは（23）に代入して、ニュートンラプソン法を適用することになるわけである。この
際の注意点は、式の形である。式（17）あるいは（23）の分母には変数を含む項がある。KM は 1 より大きい
ので、この項がかなり小さくなってしまい、上手く集束しない。分母の、変数を含む項をあらかじめ両辺にか
けて分母から消しておけば、上手く集束して、解が求まる。初期値は、あり得る範囲で取れば良い。筆者は、
初期値に 0 を用いており、それで上手く計算できている。
7
実勢の計算は、方程式を解くプログラムを利用すると良い。筆者は、
「maxima」というオープンソースのプ
ログラムを利用した*2 。その際のコードを付録（7.1 節）につけておく。
3 平衡定数の温度依存性
平衡定数は、エンタルピーや比熱などの熱力学データから計算できる。温度の何次までの関数とするかで、
式や温度定数の値が変わる。ここでは、有限会社コムテック・クウェストが提供している資料*3 を参考にした。
4 自由エネルギーを最小化する方法
自由エネルギーを最小化する方法については、松本・横川らの日本語の解説に詳しい*4 。また、CEA パッ
ケージの著者らによる解説も、参考になるだろう*5 。
以下、自由エネルギーを最小化する方法で計算する CEA パッケージについて解説する。
4.1 CEA パッケージ
CEA パッケージは、NASA の Bonnie McBride 博士と Sanford Gordon 博士のグループが開発した平衡計
算用のコードである*6 。豊富な熱力学データが用意されていて、かなりのガス反応を扱うことができる。コー
ドは、JAVA で書かれており、Windows、Mac、Unix で使用できる。Windows 版のインストールと使い方
は、東京大学准教授三好明先生の web ページ*7 に詳しい。Mac 版も同じようにインストールできるが、GUI
版の実行ファイルは上手く動かないので、ソースからコンパイルする必要がある。
使い方は、先に紹介した三好先生の解説が分かりやすいだろう。
メタン水蒸気改質反応について計算した時の、入力ファイルを参考までに載せておく。
1. problem
2.
tp
case=SC1
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”は、温度と圧力を一定にするという計算
条件を指示している。その後に、計算する温度と圧力を設定している。もしも定体積条件で計算したい場合
*2
*3
*4
*5
*6
*7
URL http://maxima.sourceforge.net/
URL http://comtecquest.com/ userdata/EqulibriumConstant&SteamReforming.pdf
T. Matsumoto, H. Yokokawa, Netsu Sokutei 19 (1992) 170―173.
S. Gordon, B.J. McBridge, Computer Program for Calculation of Complex Equilibrium Compositions and Applications: Part I. Analysis, NASA Reference Publication, 1994.
URL http://www.grc.nasa.gov/WWW/CEAWeb/
http://www.frad.t.u-tokyo.ac.jp/public/cea/cea index.html.ja
8
は、”tv”と書き、その後、温度と体積を指定するわけである。
3 行目から反応物を指定している。ここでは、ガスの名前とモル数で指定している。これは、S/C = 1 の例
である。
6 行目が大事である。ここでは、生成物として考慮するガスの種類を指定してある。これを指定しなくても
計算でき、その方がよりエネルギーの低い熱力学的安定状態を与える。しかし、2.3 説で述べたように、起こ
る反応を限定して、平衡を計算したい場合がある。その時に、この”only”を利用しよう。ここでは、メタン水
蒸気改質反応と水性シフト反応で現れるガスだけに計算対象を限定している。
8 行目は出力の形式の指定で、9 行目は終了コードである。
このように、反応条件を指定して、生成物を与え、必要に応じて反応を限定するだけで、平衡が計算できる
有り難いプログラムである。
5 比較
P0 = 1 bar
5 bar
10 bar
100
Conversion (%)
80
60
40
Present CEA
S/C = 2
S/C = 3
S/C = 4
20
0
600
図1
700
800
900 600
700
800
900 600
700
800
900
Temperature (˚C)
Temperature (˚C)
Temperature (˚C)
平衡定数から計算した平衡メタン転化率（実線、破線、一点鎖線）と CEA パッケージを利用した計
算値（丸、四角、三角）の比較
平衡定数から計算した結果と CEA で計算した結果を比較してみよう。図 1 は、種々の S/C 比、全圧 P 0
における平衡メタン転化率を計算したものである。炭素の析出はないと仮定している。また、窒素ガスな
どのバッファーガスは入れていない。平衡定数からの計算結果（実線、破線、一点鎖線）は、CEA を利用
した自由エネルギーを最小化する計算結果（丸、四角、三角）とよく一致している。ここで計算した値は、
9
Rostrup-Nielsen*8 が報じている値とも良く一致した。ただし、彼は計算の方法や前提（炭素析出やバッファー
ガスの有無）を記していないので、ここで詳しく比較の結果を述べるのはやめる。
6 参考にした文献
参考文献は、適時、脚注に載せた。それ以外の、考え方を教えてくれた、重要な 2 つの文献を挙げておく。
• Cpt. 6, ”Physical Chemistry” by Walter J. Moore.
•「化学ポテンシャルと平衡定数」小口達夫、梶本興亜、山崎勝義 著*9
7 付録
7.1 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;
*8
*9
J. Rostrup-Nielsen, in:, Phys Chem Chem Phys, 2001, pp. 283―288.
ダウンロード URL http://home.hiroshima-u.ac.jp/kyam/pages/results/monograph/Refthermo G.pdf
10
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,
11
print(round(i-273.15),float(round(conversion*10)/10))),
return("done."));
MSRvsT(600,900,10,2,1);
12