3633

チュートリアル

汎用マルチフィジックス有限要素解析ソフトウェア COMSOL Multiphysics を取り上げて新しい計算科学工学の仕事の流れを解説していきます。今回はチュートリアルの第2回目として、化学反応系を中心に解説していただきます。なお、前号よりチュートリアル記事は1ページ目のみを本誌掲載し、続きは日本計算工学会 HP 上で公開していますので、そちらも併せてご参照ください。

COMSOL Multiphysics による計算科学工学– 化学反応系（2）

橋口 真宜
1

はじめに

ていくことになり、科学技術の発展にとって大きな推進力になると思います。芸術分野や文系の方々にも計算シミュレーションに参加して欲しいものです。そうすればあたかもヒトの左脳と右脳が協働し、思いがけないヒントを得ることにもつながるでしょう。
それでは専門家以外の方が本当に計算シミュレー
ションを行えるでしょうか。答えは Yes です。
Fermi 先生のように手計算はできないかもしれません
が、私たちは身近な PC を使って Fermi 先生よりも複雑な内容のものを簡単に計算できるシステムを利用できます。統合型 GUI（Graphical User Interface）をもつ計算環境を利用すれば、問題解決のための操作も手続きもさほど複雑に感じることはありません [2]。
身近にあるノート PC やデスクトップ PC を利用しながら計算シミュレーションを行うことは計算科学工学
にとって大きな意味があります。ある研究者は、電卓でボタンを押しながら一個ずつ数字の動きを見ていったことがカオスの発見につながったのであり、もしスーパーコンピュータを使っていたら大事なことを見逃していただろうと述懐しています。
今回はチュートリアルシリーズ第2回として、化学反
応工学系の計算科学工学を見ていきます。統合型計算システムとしては COMSOL Multiphysics[3] を使っていきます。現時点での最新バージョンは5.3（2017年4月26日
日本リリース）です。

本チュートリアルシリーズはマルチフィジックス解析を通じて計算科学工学における諸現象の予測を行う
方法について解説しています。マルチフィジックス解
析は伝熱、構造、流体、化学、電磁気、量子力学、音響、生物物理、数理科学、食品科学といった広い分野にまたがるため、専門分野外の領域へも足を踏み入れる必要があり、そのような場面では自信のもてる解析をどのように組み立てていくのかが大きな課題となっています。
皆さんは Back-of-the-envelope calculation という言葉をご存知でしょうか。文字通り、紙きれなどの上で簡単な 計算 を するこ と です。 イ タリ アの 物理 学 者 Enriko
Fermi は次元解析、近似計算、推定などを行う際に、この方法が得意だったことで有名です。フェルミ問題とも呼ばれています [1]。
本 シ リ ー ズ で は、 身 近 な PC を 使 っ た Back-of-theenvelope calculation の考え方の発展形や自信のもてる解析を組み上げることも考えていきます。
マルチフィジックス解析は現象を記述することが目的ですので実験についても知っておく必要があります。本シリーズでは専門家の先生方にお願いして実験との関係についても触れるようにしていきたいと考えています。
計算科学工学は、計算シミュレーションの専門家だ
けではなく、実験に携わる方々も計算シミュレーショ
ンも行うようになれば予測と観測がより密接に関係し続きは Web で
日本計算工学会誌「計算工学（Vol.22, No.3）」HP:
http://www.jsces.org/activity/journal/

筆者紹介
はしぐち まさのり
計測エンジニアリングシステム（株）第1技術部部長。COMSOL Multiphysics の 端 緒 を 切 り 開 い た ス
ウェーデンの数学者 Germund Dahlquist 先生とともに。COMSOL 社の受付ロビーにある先生の記念碑を出張の折に撮影し、ここに一緒に掲載しました。

Vol.22, No.3 2017

注．COMSOL, COMSOL Multiphysics, COMSOL Desktopおよび Livelinkは
COMSOL ABの商標または登録商標です。その他の製品名、サービ
ス名、組織名は、各組織の商標または登録商標です。なお、本文中に™、®マークは明記しておりません。本稿の説明は著者独自の考えに基づいており他を代表するものではありません。

（25）

チュートリアル

2

COMSOL Multiphysics による計算科学工学 – 化学反応系（2）計算科学工学のイメージ

いま、数学モデルによって予測したい量を x とします。x は未知数あるいは解といいます。ここでは、取り扱う現象の数学モデルは図1に示す x の二次式で記述されるとします。皆さんは図1の二次式程度であればこの式を数値的に解くシステムとして Excel を利用できます。Goal Seek（データタブ〜What-If 分析）という方法
で 解 を 求 め る こ と が で き ま す。 実 際 の 画 面（Excel2013）での解法を図1に示します。答えは -2と出ます。
皆さんはこの二次式では因数分解ができることを知っていますので、因数分解を行い、4も解であることがわかります。Excel ではセル A1に改めて1という数値を入れた後に Goal Seek を実行しますと、解4を算出します。
アドオンでソルバーを追加すればそのソルバーを制約条件付きで実行することで4を求めることができます。

これを見れば二次式を放物線として観察でき、y=0との交点である解は x=-2と4の2個に限ることも確認できます。この図があれば、先ほどの Excel で -2とは異なる解を求めるためのセル A1の初期値は4に近い数値を入力すれば良いこともわかります。これは手計算ではありませんが、Back-of-the-envelope calculation の考え方の応用として私たちにも簡単にできるものです。計算科
学工学においては解を得るだけではなく大局的な様相を理解しておくことが大切です。
ここで示した事例では、解を求めるために Excel を利用し、図を描くのに gnuplot というソフトウェアを使いました。もちろん Excel でグラフを描くことはできます
が、一般には自分の慣れているソフトウェアを組み合わせて作業を行います。その場合、GUI（Graphical User
Interface）は各ソフトウェアによって異なります。現在
では統合 GUI 型のソフトウェアを利用すれば、マルチ
フィジックス解析全体の作業を同じルック & フィールで行えるようになっています [2]。慣れない分野のフィジックスも取り扱うマルチフィジックス解析ではここが非常に重要なポイントになります。

3

図1

x に関する二次式とエクセルでの解法

さて、x は -2と4であることは分かったのですが、これで現象が理解できたという気持ちにはならないと思います。それでは、現象を記述する図1の二次式の全体像を絵で描いてはどうでしょうか。
図1の二次式の場合にはフリーソフトウェア Gnuplotを使うと簡単に表示できます。コマンドプロンプト画
面（Windows スタートメニュー〜cmd を打ち込むと表示）
で gnuplot と入力することで、gnuplot を立ち上げます。
gnuplot > set xrange [-4:6]
gnuplot > set yrange [-10:10]
gnuplot > plot x**2-2*x-8 with lines linewidth 8と入力後に Enter キーを押せば図2のグラフが表示され
ます。

図2
計算工学

グラフによる二次式の全体像の把握

化学における計算

それでは化学における計算を見てみます。

3.1 物理化学と反応工学ある大学の先生が化学は物理化学と反応工学の二つの立場でみると良いと教えてくれました。
物理化学は、熱力学（化学反応のエネルギー論）、量子化学（分子構造）、化学反応論（化学反応の速度）の3
つの主要な領域があります。一般化学や有機化学におけるいろいろな規則は量子論から自然に導出されます。マッカーリとサイモンの物理化学（原著1997年出
版）では分子論的アプローチという副題がついている通り、分子の見地から熱力学、化学反応速度論を理解することを目指しており、計算量子化学の章も設けてGAUSSIAN94のプログラム説明まで載せています [4]。
近似や経験的なパラメータなどを含まない最も基本的な基本法則を第1原理といい、物理化学では第1原理計
算（ﬁrst-principle calculation, Ab initio calculation）を 構 築
す る こ と を 目 指 し て い ま す。 そ れ ら は ハ ー ト リ ーフォック法と密度汎関数理論に大きく分けることができます [5]。第1原理計算の実際問題への適用は身近でも行われており、例として、第22回計算工学講演会で報告されたリチウムイオン電池の劣化要因をカチオンミ
キ シ ン グ に 着 目 し た 計 算 で は Nudged Elastic Band
（NEB）法を用いた第1原理計算が活性化エネルギーの算出に使われています [6]。
一方で、反応工学の教科書を見てみますと、反応工学では反応装置の設計・操作の観点から物事を考えていくことが重要であり、反応工学は、化学反応や生物化学反応の速度過程を、物質移動、熱移動などの物理現象を考慮して解析し、その結果に基づいて反応装置を合理的に設計し、安全に操作するために必要な知識（25-2）

Vol.22, No.3 2017

チュートリアルCOMSOL Multiphysics による計算科学工学 – 化学反応系（2）を体系化した工学であるとされています [7]。反応器の設計を主眼とした立場です。
Octave Levenspiel 先 生 の 著 し たChemical Reaction
Engineering[8] の序論にも、（以下、抜粋）「Chemical reactionengineering is that engineering activity concerned with the
exploitation of chemical reactions on a commercial scale. Its
goal is the successful design and operation of chemical
reactors, and probably more than any other activity it setschemical engineering apart as a distinct branch of theengineering profession.」と明確な記述があります。
さきほどの化学の先生から Octave Levenspiel 先生に関する面白い話も伺いました。それをきっかけに調べてみますと、Octave Levenspiel 先生は紙切れの上で計算をするのが得意だったようです。このように、紙切れの上で計算をすることを Back-of-the-envelope calculationと呼びますが、Octave Levenspiel 先生の弟子である John
Eaton が数値計算システムを開発した折に、先生の名前
Octave を付けたのが有名な数値計算ソフトウェア GNU
Octave[9] だそうです。多分、思いついたらすぐに計算できるという意味を込めたのでしょう。紙の上で計算をする Back-of-the-envelope calculation はすでに述べたよ
うにイタリアの物理学者 Enriko Fermi も有名です。皆さんもレストランに行ったときなど注文品が出てくるま
でに時間がかかるようでしたら紙ナプキンの上などで問題を解いたりしたら恰好いいと思います。
反応工学では化学反応に加えて輸送現象を取り扱います。これは Bird[10] に詳細な記述があります。
化学工学での計算とはどのようなものかを概観するには代数方程式、常微分方程式による取扱いについてはすでに述べた参考図書 [7][8] に加えて、Excel を利用したもの [11]、さらに、微分方程式の取り扱いを含めたも の は Finlayson[12] が あ り、Excel, MATLAB, COMSOLMultiphysics を利用しています。
3.2 状態方程式
まずは代数方程式である状態方程式について見てみましょう。理想気体の状態方程式は皆さんも良く知っています。図3に示します。P は圧力、V は体積、n は
モル数、R は気体定数、T は絶対温度、 は比体積です。反応工学では10気圧を越える高圧下での反応も頻繁に取り扱います。理想気体の状態方程式は図3に示しましたが、分子間のポテンシャルエネルギーは考慮されていません。高圧になって分子間の距離が接近してくるとその影響を考慮する必要が生じます。有名なもの に van der Waals の 式 が あ り ま す。 そ れ を 修 正 し た
Redlich-Kwong の式を図3に示します。このように、状態方程式は図3に示すように代数方程式を構成します。
そ れ で は、COMSOL Desktop（COMSOL Multiphysicsの GUI 名）の PDE インターフェースに用意されている
ODE/DAE フィーチャーで n- ブタンの比体積を求めてみ
ましょう。

計算工学

図3

状態方程式

データとしては、
Tc=425.2[K]、pc=37.5[atm]、
R=0.08206 [liter*atm/(g*mol*K)]、
比 体 積 を 求 め る 温 度 T=393.3[K]、 圧 力 p=16.6[atm] です [12]。
COMSOL Desktop を操作する上での要点としては、
図4のように、グローバル方程式の名前に未知数である
v を入力、方程式には Redlich-Kwong の式を =0となる形
で記述する（ただし、=0は書かない）ことです。RedlichKwong の式は比体積 v を未知数とする3次式ですので、二次式の場合と同じく、適切な初期値を与える必要があります。そこで、図4のように COMSOL Desktop の関
数定義および描画機能を利用してグラフを描くと1の近傍に解があると見当をつけることができますので、それをグローバル方程式の v の初期値とします。
計算された比体積は1.5064です。理想気体の状態方程式で計算すると1.9442ですので、Redlich-Kwongは高圧下での分子間距離が小さくなっていることを示しています。

（25-3）

図4

グローバル ODE/DAE による解法

Vol.22, No.3 2017

チュートリアルCOMSOL Multiphysics による計算科学工学 – 化学反応系（2）COMSOL Multiphysics による計算科学工学 –化学反応系 (２)チュートリアル

3.3 反応速度論
化学種が反応する場合、どのような時間変化を生じかを見てみます。時間変化のみに着目するので空間0
るかを見てみます。時間変化のみに着目するので空間0次元解析になります。ビーカーの中で反応を観察する
次元解析になります。ビーカーの中で反応を観察する
イメージです。時間変化を記述する式は時間に関するイメージです。時間変化を記述する式は時間に関する
常微分方程式になります。簡単のために
A と B が順反常微分方程式になります。簡単のために
A と B が順反応速度定数 kf、逆反応速度定数 kr で可逆反応を生じ
応速度定数kf、逆反応速度定数kr
で可逆反応を生じてているとします。このとき、化学反応式、化学種 A、B
いるとします。このとき、化学反応式、化学種
A、B の
の濃度
C A 、C B の時間変化、平衡定数の関係は次のよう
濃度 CA、CB の時間変化、平衡定数の関係は次のようにになります。Aと B の濃度の和は一定であることもわなります。AとB の濃度の和は一定であることもわか
かります。
ります。
dCA
= −k f CA + krCB
dt
kr = k f /keqd(CA +CB )
=0
dt

図6

dCB
= +k f CA − krCB
dt

つまり、これは質量保存則を示しています。
つまり、これは質量保存則を示しています。
それでは平衡定数
Keq Keqを変化させたときの濃度変化をそれでは平衡定数を変化させたときの濃度変化可視化してみます。COMSOLDesktopを立ち上げ、空間を可視化してみます。COMSOLDesktopを立ち上げ、
0 次元、化学種輸送：反応工学(re)、追加、スタディ：
空間0次元、化学種輸送：反応工学（re）、追加、スタ
時間依存、
完了とします。モデルビルダで反応工学(re)ディ：時間依存、完了とします。モデルビルダで反応を右クリックし、反応を追加し、そこに
A<=>B
と書き
工 学（re）を 右 ク リ ッ ク し、 反 応 を 追 加
し、 そこにます。順方向(=>)、逆方向(<=)を合わせたという意味
A<=>B と書きます。順方向（=>）、逆方向（<=）を合わです。速度定数セクションで平衡定数指定にチェックせたという意味です。速度定数セクションで平衡定数を入れ、順反応速度定数 kf を１にします。平衡設定セ
指定にチェックを入れ、順反応速度定数kf を1にしまクションで平衡式をユーザー定義として、平衡定数にす。平衡設定セクションで平衡式をユーザー定義としKeqを入力します。グローバル定義を右クリックし、て、平衡定数にKeq を入力します。グローバル定義を
パラメタを選択後、名前に Keq、式に１を入力します。
右クリックし、パラメタを選択後、名前にKeq、式に1
初期値の設定ウィンドウで、A を 100、B を 0 にします。
を入力します。初期値の設定ウィンドウで、A
を100、
スBタ
デ ィ １ ： ス テ ッ プ １ ： 時 間 依 存 でを0にします。スタディ1：ステップ1：時間依存でrange(0,0.1,15)(単位：秒)を指定後、スタディ１を右range（0,0.1,15）（単位：秒）を指定後、スタディ1を右ク
クリックして計算を実行します。ここまでの操作の様リックして計算を実行します。ここまでの操作の様子
子を図５に示します。
を図5に示します。
得られる結果は図６(a)です。A、B の濃度は初期に急激に変化し、同じ濃度に到達します。平衡定数 Keq を0.5 にすると、図６(b)のように A の濃度は B の濃度よりも高い値に落ち着き、平衡定数 Keq を 2 にすると図

図5

反応工学における可逆変化の設定得られる結果は図6(a) です。A、B の濃度は初期に急激に変化し、同じ濃度に到達します。平衡定数 Keq を0.5にすると、図6(b) のように A の濃度は B の濃度よりも
図５ 反応工学における可逆変化の設定
高い値に落ち着き、平衡定数 Keq を2にすると図6(c) の
計算工学
計算工学

６(c)のように B の濃度が A の濃度よりも高い値になります。その理由を考えてみます。それには平衡点の解析（analysis of equilibrium points）が役に立ちます。図６(d)に示すように A と B の濃度の時間変化がと可逆変化と平衡定数の関係

ように B の濃度が A の濃度よりも高い値になります。その理由を考えてみます。それには平衡点の解析（analysis
of equilibrium points）が役に立ちます。図6(d) に示すように A と B の濃度の時間変化がともに0となる関係を調べます。すると Keq=CB/CA という関係を得ることができま
図６ 可逆変化と平衡定数の関係
す。つまり、Keq
は定常状態（平衡点）での CB と CA の比になっているためであることが分かります。平衡点のもに0 となる関係を調べます。すると
Keq =C
B /C A という解析は便利な
Back-of-the-envelope calculationです。
関係を得ることができます。つまり、Keq は定常状態（平衡点）での C B と C A の比になっているためであるこ
3.4 回分反応器と速度定数と 反応工学でバッチと呼ばれる回分操作は反応原料をが 分 か り ま す 。 平 衡 点 の 解 析 は 便 利 な
Back-of-the-envelope
calculation です。
すべて反応器に仕込んでおいてから反応を開始し、適当な時間が経過した後に反応混合物全体を取り出す操３．４
回分反応器と速度定数作法です。
反応工学でバッチと呼ばれる回分操作は反応原料をす
さて、回分反応器を用いて n- 臭化プロピル（簡単のべて反応器に仕込んでおいてから反応を開始し、適当ため化学式を A とする）とチオ硫酸イオン（化学式を Bな時間が経過した後に反応混合物全体を取り出す操作とする）の液相反応を温度310.7K 一定で行ったとしま
法です。
す。速度定数 k としたとき、CA、CB の時間変化は時間
さて、回分反応器を用いてn-臭化プロピル（簡単のたに関する常微分方程式で記述されます。反応速度rは
め化学式をA とする）とチオ硫酸イオン（化学式を-kCACB です。各時刻におけるA の濃度 CA は CB を使っBとする）の液相反応を温度
310.7 K 一定で行ったとして次式で計算できます。
ます。速度定数kとしたとき、C
A ,C B の時間変化は時間
CA=CB-57.01
に関する常微分方程式で記述されます。反応速度rは
ただし、A の初期濃度 CA0=39.58 mol/m3、B の初期濃
3－kC
B です。各時刻における度 CAB0C=96.59mol/mとします。 A の濃度 CA は CB を使って次式で計算できます。
回分器内での B の濃度は以下の通りに実測値がある
Cとします。時間は0秒から11232秒までのデータです。
A =C B －57.01
ただし、A の初期濃度 CA0 =39.58 mol/m 3 、
B の初期濃度 C B0 =96.59mol/m 3 とします。
回分器内での B の濃度は以下の通りに実測値があると-9 3 秒までのデータです。
します。時間は0 秒から 11232
この実測値からk=1.62×10
m /(mol s)と求まります [7]。

反応工学は実測値が必須であり、実験に携わる人が計算解析を行なうようになれば、そこの理解が進み、計
この実測値から
k=1.62×10 -9 m 3 /(mol s)と求まります。
算に必要な実験データも整備されやすくなるでしょう。
[7]COMSOL Multiphysics の化学反応工学モジュールの反応工学は実測値が必須であり、実験に携わる人が計
反応工学インターフェースを利用して計算する手順を
算解析を行なうようになれば、そこの理解が進み、計示します。反応工学インターフェースは、空間0次元の算に必要な実験データも整備されやすくなるでしょう。
解 析 用 に、 回 分 反 応 器、CSTR（連 続 槽 型 反 応 器、Continuous stirred tank reactor）、半回分反応器、管型反COMSOL Multiphysics の化学反応工学モジュールの反（25-4）(??-4)

Vol.22, No.3 2017Vol.22, No.3 2017

チュートリアルCOMSOL Multiphysics による計算科学工学 – 化学反応系（2）応器（プラグ流）のモデルを用意しています。ここではその中の回分反応器（バッチ）を使います。回分反応器
では良く撹拌されており空間分布を無視できるので空間は0次元と考えます。
そ れ で は、COMSOL Desktop を 立 ち 上 げ、 空 間0次元、 フ ィ ジ ッ ク ス： 化 学 種 輸 送： 反 応 工 学（re）、 追加、スタディ：時間依存、完了とするとコンポーネン
ト1がモデルビルダに組み込まれます。回分反応器はデフォルト設定になっています。反応工学（re）を右クリックして反応式を入力します。＋や＝、＞を使って記述します。その結果、入力した反応式に基づいて反応速度が ref.kf_1*re.c_A*re.c_B と自動的に算出されます。re は反応工学を、kf_1は速度定数、re.c_A、re.c_B
は CA、CB に対応しています。化学反応では左辺の反応
系の量論的係数を負としますので ν は A も B も -1です。
つまり、CA も CB もこの反応速度によって時間的に減少
していきます。続いて初期値を実測値の通りに設定します。空間0次元ですのでメッシュを作成する必要はありません。スタディ1：ステップ1：時間依存をクリッ
クして、range（0,1,11232）と入力します。単位は s（秒）です。スタディ1を右クリックして計算を実行します。
計算時間は2秒です。現象時間11232秒の計算は一瞬で終了します。
結果の下に濃度（re）が自動作成され、CA と CB の時間
変化がグラフ化されます。前述の時間と濃度 CB の実測
値をカンマで区切られた csv ファイルとして用意し、COMSOL Desktop の関数：補間機能を使って読み込みます。補間表に直接、入力することもできます。関数
名 を cB_exp と し て お き ま す。 こ れ を 使 っ て、comp1.
cB_exp(t)-57.01、comp1.cB_exp(t) を 丸 印 で プ ロ ッ ト してやると、計算結果（破線）は図7のように実測値を良く予測していることがわかります。

図7

よって発火させて燃焼させますが、HCCI では圧縮による自己発火のみで燃焼させます。ガソリンエンジンとしては究極のエンジンとも考えられマツダに見るように燃費の飛躍的向上を目指せます [13]。
メタンを利用した HCCI 燃焼反応を総括反応で表した反応式は次の通りです。
CH4+2(O2+3.76N2) → CO2+2H2O+7.52N2
実際の反応は多段階素反応によって生じます。この場合、53種類の化学種と325個の素反応式で表現されます。反応式は GRI-Mech 3.0[14] としてインターネットから入手できます。燃焼計算では素反応を使います。
これだけの数の反応式を解こうとするとシステム化が必要です。GRI-Mech 3.0は CHEMKIN のフォーマッ
トで記述されています。COMSOL Multiphysics の反応
工学インターフェースは CHEMKIN のフォーマットで反応式を読み込むことができますので多段階反応式を容易に入力できます（図8参照）。

図8

さて、燃焼反応では先ほどの液相反応とは異なって、発熱が生じ、系の温度も時間的に変化します。エンジンでは燃焼室容積もピストンの往復運動に応じて時間的に変化します。空間3次元の計算を実施すれば相当な計算量になることが予想されます。
いろいろな設計検討をするには計算時間の短い解析が適しています。
COMSOL Desktop で反応工学をモデルビルダーに読み込み、反応器タイプをバッチとして、方程式ビューを参照すると、図9の右上に記述した方程式系が表示さ

CA、CB の時間変化と測定との比較

3.5 HCCI 燃焼と空間0次元解析
HCCI（予混合圧縮自動着火、Homogeneous Charge
Compression Ignition）による燃焼反応です。通常のガ
ソリンエンジンでは圧縮した混合気を点火プラグに
計算工学

CHEMKIN データのインポート（25-5）

図9

HCCI 燃焼エンジンの0次元解析
Vol.22, No.3 2017

チュートリアルCOMSOL Multiphysics による計算科学工学 – 化学反応系（2）れます。質量保存則、熱力学第1法則、発熱量の表現式です。式中の体積 Vr をエンジン燃焼室内容積と考え、シリンダの上下動による Vr の時間変化を図9下方のように与えてやれば計算ができます。反応器容積 Vr に変数 Vol を入力し、Vol の具体的な内容は図10のように式で定義します。
数式の入力においては図10の通り、自分のノートの内容そのままの形で記述できます。図9の式と見比べてみてください。

図10

図11

エンジン燃焼室容積 Vol の式設定

化学反応は一般に温度変化を伴い、それに応じて定圧比熱、エンタルピー、エントロピーが温度依存性を持ちます。それにはまず、定圧比熱を温度の多項式で近似します。



2 a T 3 3 a T 4 4
Cp,iRg aa1 +aa22TT+a3aT3 T2 C
+4a4 T +5 a5 T
p ,i= R
aa2 2 a3a3 3 a a4 4 a5 a55 hhi= R
Rgg  aa11TT + 2 TT 2+ TT 3 + 4 T T4 + T T5 a+6 a6
22
33
44
5 5




a3 2 a4 3 a5 4
si = Rg a1 In
InT + a2 T + T + T + T + a7234

一般にはこの形をした NASA の温度多項式に従います。エンタルピー、エントロピーはその積分として算出できます。
GRI-Mech3.0には熱力学データも用意されており、
CHEMKIN フォーマットを介して COMSOL Multiphysicsに読み込めます（図8参照）。
COMSOL Multiphysics ではメタンによる HCCI エンジンの例題を用意していますのでそれを使った計算結果を示します。興味を持たれた読者は関連資料を入手できます [15]。
計 算 を 実 行 さ せ る と1.6GHz ノ ー ト PC で メ モ リ 約
1.3GB、 計 算 時 間51秒 で し た。 計 算 内 容 は、 回 転 数
1500rpm で時間0.026秒までの時間依存計算を H2O の初期モル分率を4種類変化させるパラメタスイープ付きというものです。図11に燃焼室内の圧力の時間変化の計
算結果を示します。
化学反応におけるこの0次元計算は現在の Back-of-theenvelope calculation であると考えられます。化学反応および熱力学データが GRI-Mech 3.0のように一つのセットして完備されていればこのように実用に近い場面での応用計算もすぐに実現できます。

が、実際の場ではいつもそばにいるわけではありません。従って、実際の反応器の空間形状と対で反応を考える必要があります。そのために空間依存モデルが必須です。
複数の成分気体が反応して体積膨張をする場合には混合気体の密度が変化するので混合気体は圧縮性流体として取り扱う必要がありますが、COMSOL Multiphysics ではそれもすでに用意されているものを利用できます。
簡単な例題として、化学種 A（質量分率 wA）が解離して化学種 B（質量分率 wB）になる反応を考えます。図12に反応式、必要な基礎方程式、データを示します。実
際の計算には誤差がつきものです。質量分率の和は1であることから、質量分率 wA だけを解き、wB=1-wA として算出することにします。
手順としては空間0次元で反応工学インターフェースを立ち上げ、反応で A=>2B を入力、速度定数にアレニ
ウス式を選択し、順反応頻度因子 Af に41.3、順反応活
性化エネルギーEf に30000を入力、反応エンタルピーを
ユーザー定義として反応エンタルピーに -80000を入力します。
初期値は初期温度473K、A の初期濃度を25mol/m3、B の初期濃度を22mol/m3としてスタディ1：ステップ1：
時間依存で range（0,1,10）（単位：秒）の範囲を計算しま

3.6 空間依存モデル
いままでで反応工学における質量保存則とエネル
ギー保存則の適用について見てきました。実際の反応
系は空間次元を持ちます。ビーカーや回分反応器内では化学反応を起こす相手がいつも横にいるわけです
計算工学

HCCI 燃焼圧力と初期 H2O モル分率との関係（25-6）

図12

解離反応、基礎式およびデータ
Vol.22, No.3 2017

チュートリアルCOMSOL Multiphysics による計算科学工学 – 化学反応系（2）す。うまくいけば反応部分の設定は完成です。ここまでで0次元反応速度論が完成しました。
次に、この反応を空間依存モデルに拡張します。反
応工学（re）を右クリックして、空間依存モデル生成を
モデルビルダに追加します（図13参照）。

図13

空間依存モデル作成ノードの追加

ここでは質量分率を使った計算を行うので COMSOLMultiphysics の化学反応工学モジュールの高濃度種輸送
インターフェースを利用します。空間3次元とし、高濃度種輸送、層流、伝熱（流体）、スタディタイプを定常に選択（図14(a) 参照）し、空間依存モデル生成の作成 / リ
フレッシュボタンをクリックすればコンポーネント2が生成され、指定したフィジックスと化学インターフェースが追加されます（図14(b) 参照）。化学インターフェー
スには反応工学で設定した内容が引き継がれます。各化
学種のモル質量はこのインターフェースで入力します。
3次元の円筒形流路をジオメトリで作成し、各フィジックスを同じルック & フィールの GUI 操作で必要事項を設定します。メッシュは円筒の端面に2次元の三角メッシュを生成し、円筒流路の軸方向にスイープメッ
シュを作成します。その後、円筒形側面に境界層メッ
シュを作成すれば完成です。まず、スタディ1で化学、流体と高濃度種輸送のみを選択して等温度における速度場、圧力場、質量分率の定常解を求めます。続いてスタディ2ではその定常解を使って、上記のものに伝熱（流体）を加えた計算を行い定常解を求めます。定常計
算では良い初期値（Initial guess）があればそれだけ収束性が向上します。COMSOL Desktop のモデルビルダには図15のように各スタディにおいて、解析するフィ
ジックスを選択できる機能を備えています。さらにス
タディ1で得た解を引き継いでスタディ2の計算を実行することができます。計算時間は全体で2分程度でした。
実際の適用においては実験値や既存の数値解などとの比較によって計算結果を検証する必要があります。ここで示した3次元計算結果は別途行った軸対称2次元の定常計算値と比べて妥当な解であることを確かめました。
図15に示したフィジックスの選択機能は、マルチフィ
ジックス解析を行なう上で重要な機能です。計算を構築するプロセスでは、個々のフィジックスが正しく計算されているかをまず確認し、その後、連成度を強くしていく作業を行います。この選択機能はフィジックスの
連成 / 非連成を自由に行えるので大変便利です。
計算工学(a) 反応工学からの空間依存モデルの設定(b) コンポーネント2の内容と計算結果
図14

図15

空間3次元反応流解析

複数スタディとフィジックス選択

3.7 その他の反応様式
反応工学は、いままでみてきた液相反応、気相反応のほかに表面反応が重要である場合があります。
例えば、自動車の排ガス浄化には触媒が使われてお（25-7）

Vol.22, No.3 2017

チュートリアルCOMSOL Multiphysics による計算科学工学 – 化学反応系（2）り、触媒表面上の反応機構が性能に影響を与えます。
第22回計算工学講演会では尿素 SCR（選択還元触媒、Selective Catalytic reduction）を対象とした数値計算が報告されています。そこでは多孔質媒質中の多成分の流体計算、化学種の希釈種輸送計算、温度計算の例が示されています。さらに COMSOL Multiphysics の特徴である偏微分方程式を自由に記述できる点を活かして、
触媒上でのアンモニアの吸着 / 脱離反応に関するユー
ザー独自の定式化を、解析に組み込んでいます [16]。
このような混合気の計算においては輸送特性も計算する必要があります。反応工学インターフェースには、多成分化学混合物の相互拡散係数（気体運動論による式、溶媒中の液体に関する Wilke-Chang 式）、混合物の粘性係数、熱伝導係数（Stiel-Thodos 式など）の計算式が組み込まれています。ユーザー定義も可能です。
化学工学分野では、上述の例のように、多孔質媒体が良く利用されます。
多孔質媒体内の流れを解析する場合、流速が遅い場
合にはダルシ―則が利用されます。多孔質媒体内の流速が速くなり、かつ多孔質に隣接する非多孔質領域の流れ場はナビエ - ストークス方程式によって記述される
場合には、多孔質媒体内の流れを記述する式もダルシ―則をナビエ - ストークス方程式に近づけたブリンクマ
ン方程式を利用します。
非多孔質領域および多孔質媒体領域では化学種が輸送されるため、それらを自由に取り扱える取り回しが必要です。
COMSOL Multiphysics はダルシー則を拡張したブリ
ンクマン方程式を利用するインターフェース、非多孔質・多孔質媒体中の化学種輸送を取り扱うインター
フェースを共に備えており、そのような場面に対応できます。図16はそれらのインターフェースを利用した解析例です。図16(a) に示すように、多孔質媒体中で化学種 A と B が反応して化学種 C が生成される計算です。
図16(d) を見ますと、多孔質媒体で化学種 C が生成されている様子がわかります。計算による可視化を利用す
れば多孔質媒体の位置と反応との関係も理解しやすくなります。
次回は化学反応に電子が関与する電気化学について見ていきたいと思います。

計算工学

図16

多孔質媒体のある化学反応流の計算

■参考文献
[1]
[2]

[3]
[4]
[5]
[6]

[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]

（25-8）

https://en.wikipedia.org/wiki/Fermi_problem.
橋口真宜、COMSOL Multiphysics による計算科学工学 新 し い プ ラ ッ ト フ ォ ー ム (1)、 計 算 工 学（Vol.22 No.2
2017）.
https://www.comsol.jp/products, http://www.kesco.co.jp.
マッカーリ、サイモン、物理化学（上）、東京化学同人、
1999.
白 井 光 雲、「第1原 理 計 算 と 密 度 汎 関 数 理 論」、2005、
http://www.cmp.sanken.osaka-u.ac.jp/~koun/Lecs/dft.pdf.
山中拓己、狩野恒一、高岸洋一、山上達也、3次元凝集構造リチウムイオン輸送解析による電池特性劣化評価、計算工学講演会論文集 Vol.22（2017年5月）.
橋本健治、反応工学、培風館、2009.
O C TAV E L E V E N S P I E L , C H E M I C A L R E A C T I O N
ENGINEERING, Third Edition, Wiley, 2016 (Reprint).
https://www.gnu.org/software/octave/.
Bird et al, Transport phenomena, Revised Second Edition,
Wiley, 2007.
伊東章、Excel で化学プロセス計算、丸善出版、2014.
Bruce A. Finlayson, Introduction to Chemical Engineering
Computing, Second Edition, Wiley, 2014.
https://newcars.jp/news/mazda-hcci-2018/.
http://www.me.berkeley.edu/gri_mech/.
https://www.comsol.jp/model/homogeneous-chargecompression-ignition-of-methane-1852.
佐久間優、茂木春樹、COMSOL Multiphysics® を用いた排ガス浄化シミュレーション、計算工学講演会論文集、Vol.22（2017年5月）.

Vol.22, No.3 2017