膜モジュール設計シミュレーションシステムの開発
西村 拓朗, 船津 公人
Return
1 はじめに
水素は、基幹化学物質の合成原料や石油精製、石油化学等における反応剤としての利用など、広範な産業分野において大量に使用されている。さらに、次世代クリーンエネルギー源としての水素利用エネルギーシステムに関心が集まり、実用化に対する期待が高まってきた。この社会ニーズに対して、エネルギー原単位に優れた高効率な水素製造プロセスの開発が急がれている。水素供給源としては水の電気分解によるプロセスが提案されているが、わが国では電力単価の高さが問題となっている。
本研究は下記に示す天然ガス(主成分：メタン)の水蒸気による改質反応に対応したもので、天然ガス源としては近海のメタンハイドレートが有望視されている。
CH4 + H2O ⇔ CO + 3H2 水蒸気改質反応
CO + H2O ⇔ CO2 + H2 CO変性反応
この水素製造プロセスのキーテクノロジーは原料となるメタンを水素に変換・改質する触媒反応技術と、発生した水素を不純物から分離・抽出する精製技術であり、その技術課題は、水素の生成工程と分離・精製工程の独立に起因する熱効率問題や、全体システムの複雑さによる装置の大型化である。ここで水素生成工程と水素分離工程を一体化することにより、反応温度の低減(500℃)や圧縮動力の低減が可能となり、プロセスの簡素化による小型化や、水素製造時のエネルギー効率の大幅な向上が期待されている。
本研究では、多孔質無機膜を対象とした500℃以上での化学反応プロセスを利用して水素生成と水素分離を一体的に行うことを特徴とする高効率高温水素分離膜のモジュール化技術開発に際して、分離効率の向上を目的とするガス分離膜モジュール構造最適化システムの開発を進めている。具体的には、分離膜モジュール構造の設計において、制限空間内にどのような分離エレメントをどのように配置すれば最も分離性能が上がるか、という課題を解決するのがこのシステムの目的である。
Figure 1. Hydrogen separation module.
2 システムの全体構成
ガス分離膜モジュール構造設計支援システムは、
-
分離膜の集合体であるモジュール構造データを創出する部分
-
モジュール構造と分離条件の最適化を行う部分
-
計算流体力学（Computational Fluid Dynamics: CFD）シミュレーションを実行する部分
から構成される。
Figure 2. Optimization support system in the design of the module for hydrogen separation.
当システムはWindowsXP上のVisual Fortranを用いて構築し、分離膜モジュール構造可視化ツールとしてのグラフィック表示にはOpenGLを活用した。取り扱う分離膜モジュール構造の仕様は以下の通り。
| モジュール直径 | 1,000mm以下 |
| エレメント直径 | モジュール直径以下 |
| 長さ | 2,000mm以下 |
| ピッチ | エレメント直径以上 |
| 配置できる膜エレメント本数 | 10,000本以下 |
3 膜モジュール構造データの発生
提案された構造パラメータセットにより、目的とするモジュール構造の座標データを発生し、水素回収率計算用入力ファイルを生成する。モジュール外形は円筒タイプと四角柱タイプに対応し、分離膜エレメントの配置は三角格子・四角格子に対応する等、一般の多管集合工作物（分離・抽出・熱交換器等）のエレメント配置座標が出力可能であり、CFDエンジンの入力ファイルを出力する。
3. 1 格子座標の発生
実用化されている多管集合装置では、そのエレメントは三角格子上に配置するか、メンテナンス性を重視した四角格子上に配置するのが一般的である。装置集積度の点からは三角格子上の配置が良好であるが、今回は両格子から選択して格子座標を発生させ、その格子点を中心としてエレメントを配置するように対応させた。常にモジュール断面の中心にエレメントが配置されるように格子を発生させている。
エレメントの3本配置や6本配置のように、モジュール中心にエレメントを配置させない構造座標が必要な場合は、モジュール中心から発生される単一円周上にエレメントを均一に配置させるように対応させた。
Figure 3. Arrangement on lattice and on circumference.
3. 2 エレメントの配置アルゴリズム
上記の格子間隔とエレメントの直径を比較して問題がない場合（格子間隔が充分に広い場合）にエレメントを配置する。また、モジュール隔壁との干渉を判定し、干渉のないものだけを配置させた。システム内部では、モジュール半径を1.0として定常化させ、エレメントサイズや格子間隔はモジュール半径に対する比率で管理させた。これにより、モジュールサイズやその他の構造パラメーターを変更しても、2D表示のモジュール(外周円)サイズを変更することなくエレメント配置結果を表示することが容易となる。オプションとしてモジュール壁面から一定の距離にエレメントを配置させない（空間を取る）場合にも対応させた。
3. 3 エレメント表面膜の近似
CFDによる水素の膜透過計算では、膜を平板状の平面の集合として取扱うため、エレメントの円柱形状を、直交した平面の集合として取扱う必要がある。そこで、エレメント断面の円形状を直交する直線集合で近似する場合に得られる凹凸多角形での近似検討を行った。
CFD計算の精度を上げるためには、周囲長さ（言い換えるとエレメントの膜表面積）を出来る限り円に近づけることが好ましく、かつ、面積（エレメント内部の空間体積）も円に近いことが望ましい。
円に外接する凹凸多角形近似では面積は円に比べて大きめになるが、頂点数を増すに従い円に収束する。一方、周囲は円の1.273倍で不変である。円に内接する凹凸多角形では面積は円に比べて小さめになるが、同様に頂点数の増加とともに円に収束する。この場合の外周は円より小さめなもの（四角形近似で0.9倍）から頂点数の増加に従い円の1.273倍に収束する。
これらに比べて、円周上を凹凸多角形の頂点が円周の内外に交互に分散するFigure 4のようなパターンをとると、頂点数が少ない場合でも面積は円に近似される(Table 1)。円に対してこのような凹凸多角形の座標を発生させるアルゴリズムを開発し実装した。面積は凹凸多角形の頂点数の増加にしたがって円に収束されるが、CFD計算の観点からは計算に用いるメッシュ（グリッド）の増大を招き計算効率の低下を引き起こす。計算精度と、実用的な計算速度のバランスを考慮して今回は凹凸20角形による近似を採用した。
Figure 4. Approximation of the circle by concavo-convex polygon.
Table 1. Feature of concavo-convex polygon and circle.
| concavo-convex | Area | Circumference |
|---|
| polygon | (True circle ratio) |
|---|
| 12 square shapes | 3.560 (1.132) | 8.000 (1.273) |
| 20 square shapes | 3.200 (1.019) | 8.000 (1.273) |
| 36 square shapes | 3.250 (1.035) | 8.000 (1.273) |
| 44 square shapes | 3.160 (1.006) | 8.000 (1.273) |
| True circle | 3.142 (1.000) | 6.283 (1.000) |
3. 4 GUIと構造データの可視化
構造データ発生モジュールの機能確認と、発生する膨大な構造座標データの確認を容易にするために、構造パラメーターのためのGUIと、発生データのグラフィック表示ツールを作成した。(Figure 5)
ダイアログ画面のスライダーやボタンを操作することにより、必要に応じてエレメントサイズの変更や配置間隔・配置に用いる格子の種類・モジュールサイズ等の変更を行うと、結果はリアルタイムにモジュール断面を2次元表示するウインドウや3次元立体表示（内部透視可能）ウインドウに反映され、膜エレメントの内部配置や全体の立体構造を確認することが出来る。3次元立体表示ウインドウでは、視点を自由に移動させて立体構造の外周を周回して見る（構造が自転しているように見える）ことや、構造物の内部に入って見る事、それらの自動実行も可能で、2次元表示では掴みきれない情報を容易に与える。
Figure 5. Tool for structure coordinates generating.
4 水素回収率計算モジュール
分離膜モジュールの分離性能は、単に内在する分離膜の面積に比例するわけではなく、分離エレメントの太さや長さ、隣接するエレメント間の間隔等の空間サイズに基づくものと、流入出ガスの圧力や速度に起因する拡散律速に基づくもの等から生じる、混合ガスの濃度分極の影響を大きく受ける。この問題を正確にトレースするために、汎用計算流体力学（CFD）パッケージを応用した。このCFDエンジンはモジュール構造創出部分から入力データファイルを受け取り、収束計算結果をファイル出力する。実際のモジュール構造を想定して、モジュール側壁（煙突状フランジ）からのガス流入・流出や、モジュール内部の邪魔板によるガス拡散向上と滞留等の影響も反映可能である。
CFDエンジンは、高羽[1, 2]らにより水素分離膜の透過計算用に機能拡張されたCHAM社のPHOENICSを用いた。この出力結果から水素濃度（分圧）や水素回収率を計算し、これを遺伝的アルゴリズムの評価値（適応値）として構造パラメーターの最適化に関する指標とした。
CFDを用いた計算の一例を示す。Figure 6の上図はモジュール内部に膜エレメントを9本配置した場合のモジュール内部の水素分圧分布を示している。Figure 6の中図はモジュール内部に膜エレメントを一本配置した場合の水素濃度分極の様子を示したものである。水素分離膜モジュールでは分離効率を向上させるために、モジュール内部におけるガスの流れを十分に攪拌する必要がある。Figure 6の下図はモジュール内部に邪魔板を設置し、流動状態を強制的に変化させた結果である。この計算ではガスの入り口から出口までに4枚の邪魔板を設置し、図の白矢印はガスの流体の流れを示している。このように邪魔板の設置によるガスの流れを計算し、邪魔板の配置を検討することも可能となる。
Figure 6. CFD calculation result by PHOENICS.
5 分離膜モジュール構造の最適化
構造パラメーターの最適化は遺伝的アルゴリズム（以下GA）を応用した。GAにより提案された構造パラメーターセット（入力データ）と水素回収率（出力結果）をデータベースに蓄え、次世代の水素回収率計算のための構造パラメーターセットを提案する。これを繰り返して分離膜モジュールの構造パラメータを最適化する。
GAの概略はFigure 7に示されるが、データの取り扱いや運用が規定されるわけではなく、解決すべき問題に合わせて組み合わせるため実装のレパートリーは非常に広い。以下に今回の問題解決のために用いたGAにおける染色体のコーディング、及び、各段階における処理方法の採用内容について述べる。
Figure 7. Flow of a genetic algorithm.
5. 1 染色体の設計
遺伝的アルゴリズム（GA）では、一つの個体（構造パラメーターのセット）を一つの染色体で表す。言い換えれば個体の特性は染色体を用いて決定されるので、アルゴリズムを実施する前に染色体を問題に合うように設計する必要がある。
今回は構造パラメーターセットを表現するためにTable 2に示す染色体[3]を設計した。取り上げた構造パラメーターは、モジュール長さ・モジュール直径・膜エレメントの直径・膜エレメントの配置間隔（ピッチ）である。染色体の長さは、各パラメーターの探索範囲のステップ数（ｍｍ単位）を表現するのに必要なビット数を順次並記したものである。これは構造パラメーターの表現に関して不必要なビットを含まず、後述の遺伝的操作の効率を悪くしたり無意味となることを排除した結果である。
Table 2. Chromosome expression of one individual.
5. 2 グレーコーディング
GAでは問題を効率よく処理するには染色体のコーディングが重要な要素となる。バイナリコーディング（2進表記）では10進表記で隣り合っている値のハミング距離（遺伝子の異なる部分の数）が1にならない、つまり似ていない場合が頻繁に存在する。
例えば10進表記の15は2進表記では01111となるのに対して十進表記の16は二進表記では10000であり、ハミング距離は5となる。(Table 3)
このような場合、16が最高値で15がその次に評価が高いとすると、GAで生まれた子供は親に比較的に似るので、ほどほどに評価の高い15や14の周辺に評価の集合が集まっても、16は15との近似度が低い（ハミング距離が大きい）ので見つけにい。
この問題を解消するために、グレーコーディングを採用した。グレーコーディングを行うと10進数で隣り合う数字のハミング距離は必ず１になる。
Table 3. Correspondence of binary code and gray code.
| decimal | binary code | gray code | decimal | binary code | gray code |
|---|
| 0 : | 000000 : | 000000 | 16 : | 010000 : | 011000 |
| 1 : | 000001 : | 000001 | 17 : | 010001 : | 011001 |
| 2 : | 000010 : | 000011 | 18 : | 010010 : | 011011 |
| 3 : | 000011 : | 000010 | 19 : | 010011 : | 011010 |
| 4 : | 000100 : | 000110 | 20 : | 010100 : | 011110 |
| 5 : | 000101 : | 000111 | 21 : | 010101 : | 011111 |
| 6 : | 000110 : | 000101 | 22 : | 010110 : | 011101 |
| 7 : | 000111 : | 000100 | 23 : | 010111 : | 011100 |
| 8 : | 001000 : | 001100 | 24 : | 011000 : | 010100 |
| 9 : | 001001 : | 001101 | 25 : | 011001 : | 010101 |
| 10 : | 001010 : | 001111 | 26 : | 011010 : | 010111 |
| 11 : | 001011 : | 001110 | 27 : | 011011 : | 010110 |
| 12 : | 001100 : | 001010 | 28 : | 011100 : | 010010 |
| 13 : | 001101 : | 001011 | 29 : | 011101 : | 010011 |
| 14 : | 001110 : | 001001 | 30 : | 011110 : | 010001 |
| 15 : | 001111 : | 001000 | 31 : | 011111 : | 010000 |
グレーコードはnビットの２進コードが与えられたとき以下の手順で求めることが出来る。
-
2進コードの最上位ビットを、そのままグレーコードの最上位ビットとする
-
グレーコードの n-1ビット目以降を、２進コードのiビット目と i-1ビット目との排他的論理和(Exclusive OR：XOR)を用いて決める。
(i=n, n-1, ・・・, 2)
5. 3 遺伝子の選択
次世代に残す染色体を選別するのにルーレットルールを採用した。これはFigure 8に示すように適合度に比例した面積を有するルーレットをつくり、それを回して矢の当たった場所の個体を選択するという方法である。この結果、適応度に比例した割合で個体を選択することができる。つまり、適合度の高いものは高い確率で生き残るし、適応度の低い個体でもそれなりの確率で生き残れるようになる。よって、適合度に基づく淘汰・増殖を行いながら、個体集団中の多様性を維持することができる仕組みとなっている。
Figure 8. Roulette rule.
5. 4 エリート戦略
ルーレットルールでは確率的に次世代の母集団を決定するため、最良の個体が次世代に残らないことがある。これを避けるための手段として、エリート戦略を採用した。エリート戦略とは、個体の集団の中で最も適合度の高い個体（エリート）を無条件に次世代に残すという方法で、これにより最良の個体が次世代に保持され、良い個体が増加することにより局所探索能力が向上する。
5. 5 交叉
選択により適合度の高い個体が母集団内に増えるが、それだけで進化が進むわけではない。交叉により、個体同士の染色体の一部を組み替えることで、より最適な解が得られる可能性が高くなる。母集団のうち交叉に参加する個体の割合は交叉確率によって定めた。
交叉の方法は、一点交叉を採用した。一点交叉とは、染色体上のランダムな位置一つに交叉点を定め，交叉点を境に２個体の染色体情報を交換する手法である。
Figure 9. One-point crossover.
5. 6 突然変異
交叉により生ずる個体には、必ず親の個体の一部が存在するため良好な解が得られず、局所解に陥る場合がある。そこで、局所解から積極的に脱出させることを目的として突然変異操作を実施する。突然変異は,染色体上の全てのビットを確率的に反転させることによって発生させ、その結果、母集団内に新しい遺伝子を持つ個体が発生する。各遺伝子が突然変異を受ける割合を突然変異確率と言い、1/染色体長が一般的であり、これは確率的に染色体一本につき一箇所突然変異が現れることを意味する。作成したGAエンジンの操作ダイアログ画面(Figure 11)から、突然変異確率は0から1の範囲で自由に設定でき、具体的には各個体の染色体を表す全てのビットに対して乱数を発生させ、その値が突然変異確率より小さい場合にビット反転を行わせた。これにより染色体全領域での均一な突然変異を保障した。
5. 7 データベース
GAによる最適化において、適応値を算出するためのCFDエンジンの計算負荷が大きい（1セットの構造パラメータにつき数時間必要）ために、システム全体のパフォーマンス低下が懸念される。そこで、過去に計算の行われた染色体（構造パラメーターセット）を管理し、GAの中で提案された個体が既に計算済みである場合はCFDの処理に入らずに、管理されているデータベースから計算結果を呼び出す処理を加えた。GAによる最適解探索においては、極値付近で親と同様の染色体を提案する場合が増すので、適応値算出に関する計算負荷を大幅に軽減することができる。
さらにこのデータベースを利用し、複数のCFD計算機の援助を得て、CFD計算負荷を分散させることを可能とした。
Figure 10. Distribution model for Calculation load of Fitness.
5. 8 ＧＵＩとＧＡ評価結果のグラフ化
水素分離膜モジュール構造の最適化を実現するに当たって、各構造パラメーターの探索範囲の設定や評価項目の選択、及び、GAに関するパラメーターを管理するGUI（Figure 11）を作製し、ここから今回作製した構造パラメータ発生ツールやGAエンジン・CFDシミュレータを随時呼び出す統合環境を構築した。このダイアログには、最大評価結果を与えた構造パラメータを提供する部分と、GAにより評価を終えた全データを閲覧できるリストボックスを設けた。
Figure 11. Dialog for structure optimization.
構造最適化を実施中の膜モジュール構造は常に３D表示窓（Figure 12左）により確認でき、GAによる適応度と到達世代の関係はグラフ化されて別窓（Figure 12右）に表示される。これらにより、生データを閲覧することなく最適化の状況を視覚的に確認できる。
Figure 12. Result by generation advance.
6 結果と考察
今回開発した構造最適化エンジン（遺伝的アルゴリズムの応用）の性能確認のため、CFDによる水素回収率の代わりに分離膜の表面積を適応値として利用し、各条件においてTable 4の探索範囲で500世代（20個体／世代）まで実行させ、最適化のための運用検討を行った。実行環境は2GBのメモリを実装したXeon (2.8GHz) 2CPUマシンである。
Table 4. Search range of structure optimization.
| Structure parameters | Range(mm) | Step |
|---|
| Diameter of module | 50~300 | 250 |
| Length of module | 300~1000 | 700 |
| Diameter of element | 3(Fixed) | 1 |
| Pitch of element | 4~20 | 16 |
Crossover probability: 0.25（One-point crossover)
Mutation probability: 0.01
Gray code expression adoption, Roulette rule adoption
6. 1 エリート戦略の効果
遺伝的アルゴリズムの処理の中で、次世代のパラメータセットを決定する際に、前世代の最高適応値を示したパラメータセットを必ず引き継がせるというエリート戦略を用いた場合の効果を確認した。（Table 5）エリート戦略をとらずに交叉と突然変異のみに頼った場合は最適化効率が悪く、500世代まで実行しても理論最大適応値に至らない場合があることから、当問題解決にはエリート戦略の採用が必須と判断した。
Table 5. Effect of elite strategy.
| Run | Attainment generation to the theoretical maximum fitness (Individual number) |
|---|
| without elite strategy | with elite strategy |
|---|
| 1 | 369( 9) | 74(12) |
| 2 | 283(16) | 262(13) |
| 3 | 464(11) | 166(13) |
| 4 | 427( 7) | 45( 9) |
| 5 | 497(10) | 106(8) |
| av. | 408 | 131 |
6. 2 初期データセットの影響
Table 5において、実行毎に到達速度が異なるのは、第一世代の構造パラメータセットの与え方（乱数により発生）に基づくことが判明した。よって、作為的操作を排除した上で妥当な初期データセットを提案するアルゴリズムを構築した。
6. 3 計算量負荷の軽減効果
Table 4の範囲で構造パラメータの全空間探索を実行すると280万回のCFD計算が必要で、CFD計算一回あたり１時間費やすとして320年の連続計算量に相当する。GAの有効利用により、50世代（20個体/世代）の計算で最適化が達成できるとすると計算負荷は約1/3000（40日間の計算量）となり、大幅な効率化が達成できる。
さらに、GA運用に関してデータベース管理機能を活用して評価済みの構造パラメータセットのCFD計算をスキップする場合、CFD計算量を1/8に軽減できることを確認した。(Table 6)
加えて、複数のCFD計算機によりCFD計算負荷を分散させる場合は、追加計算機の有効台数が世代あたりの個体数に比例し、20個体/世代の場合で、さらに1/3から1/5の負荷削減が見込める。
Table 6. Curtailment of Fitness calculation by GA Database.
| Run | A | B | C | B/C | A' | B' | C' | B'/C' |
|---|
| 1 | 74(12) | 1472 | 207 | 7.11 | 124 | 2480 | 275 | 9.02 |
| 2 | 262(13) | 5233 | 625 | 8.37 | 312 | 6240 | 709 | 8.80 |
| 3 | 166(13) | 3313 | 415 | 7.98 | 216 | 4320 | 498 | 8.67 |
| 4 | 45 (9) | 889 | 182 | 4.88 | 95 | 1900 | 273 | 6.96 |
| 5 | 106 (8) | 2108 | 261 | 8.08 | 156 | 3120 | 360 | 8.67 |
A: Attainment generation to the theoretical maximum fitness (individual number)
A': A+50 generation
B: Number of times of fitness calculation to A (with duplication)
B': It's the same as the left. (to A')
C: Number of times of fitness calculation to A (with no duplication)
C': It's the same as the left. (to A' )
7 まとめ
遺伝的アルゴリズムによる分離膜モジュール構造最適化のためのエンジンが、モジュール形状やエレメントのサイズや配置に関して効率よく最適化できることを確認した。今後、水素分離膜モジュールの運転効率最適化における関連パラメータの包括的な最適化を進めるとともに、遺伝的アルゴリズムにおける運用の柔軟性を活用して、更なる効率化を追求する。
本研究は、「（革新的温暖化対策技術プログラム）高効率高温水素分離膜の開発プロジェクト」の一環として、NEDOから委託を受けて実施したものであり、関係各位に深謝する。
CFDエンジンの提供を受け、当システムへの組み込みに際してCFD全般について議論と助言をいただいた東京大学の中尾真一教授と高羽洋充助手へも謝意を表す。
参考文献
[ 1] H. Takaba, T. Nishimura, K. Funatsu, S. Nakao, Trans. Mater. Res. Soc. Jpn., 29, 3279-3282 (2004).
[ 2] 高羽洋充, 中尾真一, ガス分離膜モジュール設計用CFD計算モジュールの開発, 化学工学会年会研究発表講演要旨集, 640 (2003).
[ 3] Z. Michalewiz, Genetic Algorithms + Data Struc- tures = Evolution Programs, Springer-Verlag (1992).
Return