Panasonic Technical Journal Vol. 63 No. 1 May 2017

燃料電池マルチスケールシミュレーション
Multiscale Simulation for Fuel Cell

山本恵一

Keiichi Yamamoto

要藤田悠介

Yusuke Fujita

安本栄一

菅Eiichi Yasumoto

原靖

Yasushi Sugawara旨
燃料電池の理想的な触媒層構造を設計することを目的に，触媒層を構成するカーボン担体とその表面にある触媒粒子をアイオノマーで被覆した3相界面構造に基づく，商品電極サイズの燃料電池の発電性能予測できるマルチスケールシミュレーション技術の開発を行った．今回，分子動力学法によるアイオノマー内のナノスケールの物質輸送特性の算出や，カーボンの担体粒子サイズを最小単位とする多孔質構造の実効的な自己および相互拡散係数とプロトン伝導度の決定をメソスケールのシミュレーションにより行い，触媒層における電気化学反応を適切にモデル化することで，発電電圧の予測値として，実験値との乖離（かいり）±10 mV未満の精度を実現した．
Abstract
To design the ideal structure of a catalyst layer of a fuel cell, we developed multiscale simulation techniques that enable us tosimulate the performance of a fuel cell having an electrode area of a commercial fuel cell system. This fuel cell is based on a structureof a three-phase interface in a catalyst layer comprised of a carbon support and particles of catalyst covered by an ionomer. So far, wehave derived the nanoscale transport properties using a molecular dynamics method and calculated the effective self and mutualdiffusivities using the mesoscale simulation method. We consequently improved the accuracy for each scale separately. This time, weengaged in integrating a nanoscale to macroscale simulation by properly modeling the input and output data. As a result, we achieved
an accuracy in which the difference between experimental and simulation results of cell voltage was less than 10 mV.

1．はじめにスタック

固体高分子形燃料電池（以下，燃料電池）は，水素とセパレータ

モジュール
流路

空気中の酸素を電気化学的に反応させ，電気と熱を作る

MEA

クリーンなエネルギーシステムであり，低出力でも発電効率が高いという優れた特長をもち，家庭用燃料電池エカーボンの担体粒子

Pt 触媒 アイオノマー

ネファームや燃料電池自動車などに採用されている．しかしながらその本格普及には，コストが依然課題となっH次に燃料電池の構成と作動原理を簡単に説明し，開発課題について述べる．第1図に示すように，燃料電池はAir（O2）

H2ており，さらなる技術革新が必要とされている．

+H+MPL

MPL

基材

セパレータとMEA（Membrane Electrode Assembly）から構成されるセルの基本単位がモジュールとしてスタックの中で積層されている．さらにMEAは両極の触媒層およびMPL（Micro Porous Layer）と基材からなる拡散層と中触媒層

拡散層

第1図

電解質膜

触媒層

拡散層

燃料電池の構成

Fig. 1 Construction of fuel cells

央の電解質膜で構成され，各極には水素と酸素が供給され，水の電気分解の逆反応により発電する仕組みとなっなっている．ている．この発電に直接寄与する箇所が触媒層であり，そこで今回，カーボンの担体粒子サイズを最小単位と厚み約10 μmの多孔質構造からなり，粒径が数nmのPt触する触媒層の理想的な多孔質構造を設計することで，触媒，Pt触媒を担持し電子を伝導する直径数10 nmのカーボ媒層内の物質輸送特性を改善して触媒層の高活性化を図

ン担体，Pt触媒までプロトンを伝導するアイオノマー，ることを目的に，触媒層の各種構造に基づく商品電極サ水素や酸素などのガスを拡散させる細孔で構成されていイズの燃料電池の発電性能を高精度に予測するマルチス

る．これら燃料電池の構成のなかで，特に触媒層のPt触ケールシミュレーション技術の開発を行った．本シミュ媒が高価であることが低コスト化への障害となっている

レーションはナノ，メソ，マクロスケールから構成され，ことから，使用量を減らしながらも，発電効率を高める

それぞれ1 nm〜10 nm，10 nm〜10 μm，10 μm〜のオーダこと，つまりは触媒層の高活性化が筆者らの開発課題とーとなっている．この枠組みのなかで，ナノスケールに72

73

先端研究特集：燃料電池マルチスケールシミュレーション

おいては分子動力学法（MD：Molecular Dynamics Method）

分子数＋ヒドロニウムイオン数）/スルホン酸基数」と定によるアイオノマーの物質輸送特性を特徴付ける拡散係義した．過去の検討ではNafionに生体高分子の汎用ポテ数，溶解度およびプロトン伝導度の算出を行い，メソスンシャルとして有名なGAFFを使用していたが，密度がケ ー ル に お い て は 格 子 ボ ル ツ マ ン 法 （LBM：Lattice

実験値を再現しないという問題があった．それを解決すBoltzmann Method）を多成分ガス拡散対応の格子ボルツるため，本計算ではNafionにDREIDING力場をベースとマ ン 法 （MC-LBM：Multicomponent Lattice Boltzmannしたモデル[1]を採用した．さらに，高分子のねじれを表Method）に拡張するなどして，カーボンの担体粒子サイすポテンシャルでは，ファンデルワールス相互作用・クズを最小単位とする多孔質構造の実効的な自己および相ーロン相互作用に由来する非結合ポテンシャルを0に設互拡散係数とプロトン伝導度の算出に取り組んだ．また，定することで相互作用を弱めている．ヒドロニウムイオ

マクロスケールにおいては，触媒層内の物質輸送特性のンにはclassical hydroniumモデル，水分子にはSPC/FWモ改善と発電効率の向上を関連付けることを目的に，ナノデル，酸素分子にはDutkiewicz，Banerjeeのモデルを使用

およびメソスケールの結果を適切に利用することで，商した．酸素-溶媒間にはOhmoriらのポテンシャル[2]を設

品電極サイズの発電性能予測ができるようにした．以下，定し，他の異種原子間相互作用はローレンツ−ベルテロ

スケールごとに取り組んだ内容を説明する．

則により決定した．また，3次元周期境界条件下のセルにNafion・水分子・ヒドロニウムイオンをランダムに配置し，温度・圧力を与えることで含水率λ＝3，6，11，16

2．ナノスケールシミュレーション

における，実験結果と比較できるバルクNafion構造を作2.1 概要成した．高分子を取り扱うMD計算では高分子の鎖が複ナノスケールのシミュレーションにはMD法を用いた．

雑に絡み合い，系が安定構造に至りにくい．そのため，MD法では粒子間に相互作用を設定し，個々の粒子の運温度の上昇・下降を繰り返すアニーリングにより絡み合動方程式を数値的に解くことで，経時的な粒子の挙動をいを解き，T=300 K，P=0.1 MPaにおける緩和構造を得た．

（注1）

直接追うことができる．MD計算にはJ-OCTAおよびLAMMPS（SANDIA国立研究所）を用い，全原子シミュ2.3 検証レーションを行った．

まず，含水率に対するバルクNafionの密度を算出した．

今回，筆者らは最も代 表的 なアイオ ノマ ーである
Nafion（注2）

GAFFを使用した従来の計算に比べ各含水率における密を計算対象とし，シミュレーションで得られ度が上昇し，実験値との誤差3 %以内を実現した．新たるバルクNafionへの酸素溶解度を実験値と比較することに採用したポテンシャルの効果および, ねじれ角におけで，シミュレーションモデルの妥当性を検証し，アイオる非結合相互作用をなくすことにより，高分子の鎖がノマー内の物質輸送特性の算出ができるようにした．

GAFF使用時に比べ柔らかくなったため，密度が上昇したと考えられる．
次に，系に分子を仮想的に挿入し，そのエネルギー変

2.2 バルクNafionモデル
Nafionはフッ素樹脂（PTFE）の主鎖にスルホン酸基の化から溶解度を算出するテスト粒子挿入法により酸素溶側鎖をもつポリマーである．本計算で使用したNafionの解度を算出した．含水率に対するバルクNafionにおける分子構造を第2図に示す．

酸素溶解度を第3図に示す．算出された酸素溶解度は含│−OCF2ーCFCF3−O−CF2CF2−SO3 H
第2図+

Nafionの分子構造

Fig. 2 Molecular structure of Nafion

Nafionは含水することでその形を変えるため，含水量

酸素溶解度(at STP) [cm3]1.E-05−[[CF2CF2]7−CFCF2]10−

1.E-06[Sim.] Panasonic 300 K[Exp.] Gode 298 K[Exp.] Gode 333 K

1.E-07[Exp.] Parthasarathy 298 K[Exp.] Buchi 298 K[Exp.] Jeffrey 308 K

1.E-08もパラメータとなる．含水の度合いは「含水率λ=（水0510
含水率 [%]

1520（注1）（株）JSOLの日本および米国における登録商標．

第3図（注2）The Chemours Company FC, LLC.の登録商標．

Fig. 3 Water-content dependence of oxygen solubility in bulk Nafion

バルクNafionにおける酸素溶解度の含水率依存73Panasonic Technical Journal Vol. 63 No. 1 May 201774

水に対して文献[3]-[6]で示された実験値と同様に減少傾

空間補間することで，単成分から多成分系に拡張するこ向を示し，さらにこれらのばらつき内に収まることを確とを試みている．しかしながらMC-LBMの実際の計算で認した．

は対象により，数万〜数百万回程度の反復を繰り返すことから，この空間補間の繰り返しで生じる数値誤差は実3. メソスケールシミュレーション際には許容できなくなる．そこで筆者らは，並進後の空間補間は採用せず，各粒子が各ステップで確実に格子点3.1 概要に乗るよう粒子ごとに時間間隔を設定すること，およびメソスケールのシミュレーションには，内製の異種粒子間の衝突時に相手の分布関数を時間的に内挿補

MC-LBMを用いた．まず，LBMとは流体の構成粒子に対間して異種粒子間の分布関数を逐次情報交換させることする離散化された速度の分布関数を逐次的に計算し，得で多成分ガスへの拡張を実施した．

られた分布関数に対し平均化などの統計処理を行い，速度などのマクロな物理量を決定する計算アルゴリズムで3.3 出口境界アルゴリズムある[7]．LBMには，衝突項のモデルが幾つか提案されて2つ目は出口境界に関するアルゴリズムである．ここでいるが，最もシンプルで，かつナヴィエストークス方程

出口とは，ガスの平均的な流れが3次元解析領域から流出式の解を正確に再現することが知られている格子BGKする面を意味する．単成分ガスのLBMでは出口の圧力を（Bathnagar-Gross-Krook）法が現在主流となっている．

決めれば出口の分布関数は一意に決まるが，多成分ガスシミュレーション対象がガスの場合，ナヴィエストーク

ではモル分率の自由度が残るために分布関数が一意に決ス方程式が6 μm程度で適用限界に達するのに対し，定されないという問題が生じる．この問題の最も簡単なLBMは平均自由行程以下の壁との衝突も考慮できるこ

解決方法は，圧力を固定したままモル分率だけをシミュとから数nmのオーダーまで有効であると考えられておレーション領域から出口へ外挿補間することであるが，り，細孔径が100 nmを下回るような燃料電池の拡散層や実際にそれを試みた結果，この方法では数値的に安定な触媒層の移流拡散をシミュレーションするような場合に解を得ることはできなかった．これは出口部の解析領域最も有効なアルゴリズムと考えられている．しかしなが外の分布関数の情報が，解析領域内に伝達することによら筆者らの開発ターゲットである燃料電池では，各極でるものと推定される．そこで筆者らは実際の物理系にお複数成分のガスが存在することから，単成分ガスのLBMいて異種ガスが混合した場合に，短時間で各ガスの平均を多成分ガスに拡張したMC-LBMのアルゴリズム開発速度が同一化するという物理的な特性に着目し，出口にを行う必要があった．

おける全ての異種ガスの速度の法線成分が常に等しいと

MC-LBMのアルゴリズムは並進と衝突からなり，異種粒子同士の衝突の存在が最も特徴的な部分である．異種いう条件を課すことを見いだし，分布関数を出口で一意に決定することに成功した．

粒子の衝突アルゴリズムは2003年のLuoら[8]により報告されているが，筆者らは1956年のGrossら[9]の物理的仮定

3.4 検証が明確な古い分子モデルの論文を参考にしつつ，マクロ

MC-LBMのアルゴリズムとコードの検証は，理論面とスケールでマクスウェル・ステファン方程式に収束する実験面の両方から行った．理論面では，市販の汎用熱流よう異粒子の衝突演算式を確定した．3次元MC-LBMの体ソフトウェアFluent（注3）を利用して，MC-LBMの計算開発過程において，学術論文や特許などでは確認されて結果と比較し，実験面では複数の多孔質材料に対してガいない多成分ガス環境における2つの新アルゴリズムをー レ ー 透 気 試 験 （ 注4） を 行 い ， 浸 透 係 数 を 算 出 し て開発したので，次にその説明を行う．

MC-LBMの計算結果と比較した．これらの検証により理論，実験の両面からMC-LBMの十分な計算精度が確認さ3.2 並進運動アルゴリズムれたので，以下にその説明を行う．

1つ目は並進に関するアルゴリズムである．複数成分のMC-LBMはメソからマクロスケールに適用できるアガスが同時に存在する場合，軽い粒子が格子点から格子ルゴリズムであるが，Fluentはマクロスケールの汎用流点に移る時間間隔では重い粒子は格子点まで移動できな体ソフトであるため，両者の比較が可能なように，第4いことより，単成分ガスのLBMのアルゴリズムは，そのままでは多成分ガスには適用できないという問題が生じ（注3）ANSYS, Inc.の登録商標．

る．それに対して，例えば2007年のJoshiら[10]の論文で（注4）対象材料に規定圧力を掛け，規定量の空気が抜けきるは，並進後の重い方の粒子の分布関数をバイリニア法で74までの時間を測定する試験．

75

先端研究特集：燃料電池マルチスケールシミュレーション

図の右上に示すようなミリメートルオーダーの流路を模良好な一致を確認できた．

擬した構造を用いて検証を行った．このときのガスの流速は窒素，酸素ともに0.1 m/sとした．第4図の中央にモル3.5 実効的な物質輸送パラメータの算出

分率の切り出しラインと計算結果を示すが，モル分率に

次に筆者らは，MC-LBMなどを用いて燃料電池の性能ついて両者のグラフがほぼ完全に重なることを確認できを正確に予測するために必要な燃料電池の各部材の実効た．また別途，MC-LBMの結果がマクスウェル・ステフ的な物質輸送パラメータを算出した．このパラメータにァン方程式をどの程度満たすかを確認したとろ，99.5 %は，多孔質材料の浸透係数または実効的な自己拡散係数，以上の精度で同方程式を満たすことを確認した．

実効的な電気伝導度，実効的な熱伝導度，実効的な相互拡散係数などが含まれるが，本稿では燃料電池の性能への影響が特に大きいと考えられる多孔質材料の実効的な1 mm

相互拡散係数と触媒層のプロトン伝導度の算出方法を簡

モル分率 [-]

1.0

単に示す．

0.8

実効的な相互拡散係数を算出するため，筆者らは最初0.6に，多孔質材料を2本の流路で挟み，その両サイドから異0.4なる種類のガスを流入させる構成でMC-LBMによるシFluent× MC-LBM

0.2

ミュレーションを実施し，流路の出口部分の各ガスのモ0.0050

100150200

ル分率を実験的に決定することを試みたならば，燃料電

メッシュ番号 [-]第4図ル分率を算出した．もしここで，出口部分の各ガスのモ池の触媒層やMPLなどでは基材に塗工することで形状FluentとMC-LBMの計算結果の比較を保持することから，大部分の部材で実験系を開発する

Fig. 4 Comparison of the simulation results of fluent

ことができなかったと予測される．よって，この点からand MC-LBM

も，燃料電池のシミュレーションにおけるMC-LBMの有次に多孔質材料の浸透係数に関して実験とシミュレーション結果を比較した．ここで浸透係数とは，ダルシー（注5）効性が確認できる．
次に，Fluentを用いて，多孔質構造の部分に均一な相における比例定数で多孔質材料の透過の程度を決互拡散係数を設定し，それをパラメータとして変化させめるパラメータである．マクロな流体シミュレーションながら，出口部分の各ガスのモル分率がMC-LBMの結果

では浸透係数は物性値として入力するものであるが，と等しくなる相互拡散係数の値を算出した．そして，そMC-LBMでは多孔質材料の下端部分からガスを流入されを検討部材の実効的な相互拡散係数と定義した．こうせ，下端から上端にかけての圧力損失を計算することでして算出した相互拡散係数は，多孔質材料が存在しない浸透係数の値を算出できる．ここでは空隙率が異なる拡バルクの0.1〜0.8倍程度となることが分かった．一般に，

則散層のメディアを4種類用意し，ガスには空気を用い，流燃料電池のシミュレーションの分野においては，実効的速をさまざまに変化させて浸透係数の値を算出した．まな拡散係数は空隙率に比例するモデルが採用され，バルた，それに対応するガーレー透気試験を実施して各多孔クの0.6〜0.8倍程度の値がよく用いられるが，今回の筆者質材料の浸透係数の値を見積もった．第1表にMC-LBM

らの算出結果は，空隙率に単純に比例するだけの実効的とガーレー透気試験の結果を示すが，全ケースで両者のな相互拡散係数のモデルは拡散性を実際よりも高く見積もる可能性があり，十分な精度を有しないことを示唆し

第1表

浸透係数の計算と実験結果の比較た．

Table 1 Comparisons of permeability in simulation and

さらに触媒層内の実効的なプロトン伝導度を算出した．

experimental results

第5図に示すように触媒層の模擬構造に一様にアイオノ
材料A

材料B

材料C材料D

試験結果 [m2]1.52 E-15

1.15 E-15

6.8 E-12

5.7 E-12

解析結果 [m2]1.49 E-15

1.17 E-156.6 E-12

5.6 E-120.98

1.010.970.98

精度（解析/試験）

マーを被覆させ，そこに局所的なプロトン伝導度を設定し，多孔質材料の上下の境界面に適当な電位差を設定し，多孔質材料を流れるプロトンのフラックスを算出した．
そして，このフラックスを電位差で除することで実効的なプロトン伝導度と定義した．本算出の結果，局所的な（注5）圧力と，それに応じて生じる流体の平均流速を比例関係で記述する法則．

プロトン伝導度の1/15程度となる結果が得られ，具体的には1.0 S/m程度で，交流インピーダンス法などの実験か75Panasonic Technical Journal Vol. 63 No. 1 May 201776

ら推定される値とほぼ等しい結果になることが確認されとを確認した．また，定量的な検証は現在も進行中であた．これまでもプロトンの物質輸送は，薄膜のアイオノるが，電流-電圧特性に加え，酸素利用率特性，圧力損失

マー中を伝導すると想像されてはきたが，あらためて今特性，I/C依存特性，温度分布，過電圧分離などの別の評回，その定量的な裏づけを得るに至った．

価においても，高い精度で実験を再現する見込みである．

()

 αF

Vopen −circuit − φe − + φH +  （1）∇ 2 φH + = AC O2 exp  RTC
 αFODHeff+ ∇ 2 CO2 = A 2 exp Vopen−circuit − φe− + φH + 4F
 RT


()

･･････････（2）σ Heff＋：プロトン伝導度，DHeff＋：拡散係数
φe− ：電子電位，φH +：プロトン電位，

CO2：酸素濃度，A, α, F, T, Vopen-circuit：定数
第5図

走査型電子顕微鏡画像とアイオノマーを被覆した

1.1

模擬構造Fig. 5 Scanning electron microscope image and catalyst layer
電圧[V]covered by ionomer

以上の取り組みにより，実効的な自己および相互拡散1.0Exp.

0.9

Sim.

0.80.7係数とプロトン伝導度といった多孔質材料の性能を特徴0.6

付ける重要な物質輸送パラメータの算出方法の妥当性を0.50.4確認できた．

0.04. マクロスケールシミュレーション

第6図

0.20.4電流密度[A/cm2]

0.60.8

実験とシミュレーションの結果の比較

Fig. 6 Comparison of experimental and simulation results

マクロスケールのシミュレーションでは，商品電極サイズ（A4サイズ程度）の燃料電池の発電性能予測を行った．燃料電池の触媒層とMPLは50 nm〜100 nm程度の細5．まとめ孔をもつ多孔質材料であることから，発電性能予測のシミュレーションにおいては，ナノおよびメソスケールの今回，ナノ，メソ，マクロスケールの物理法則が異なシミュレーションで算出された燃料電池材料の実効的なる階層ごとに燃料電池のシミュレーション技術を構築し，物質輸送パラメータを利用する必要があった．また，A4非常に高い精度で発電性能を予測することに成功した．

サイズの電極の場合，流路の入口から出口にかけガス濃

今後は本シミュレーションを利用し，まずは物質輸送特度が大きく変化することから，触媒層とMPLの界面の酸性に注目して，カーボンの担体粒子サイズを最小単位と素濃度を入力として，触媒層における電圧のドロップ量する触媒層の多孔質構造を最適化する予定である．そのを出力する反応モデルを構築する必要があった．

後は，先述のナノスケールのシミュレーション技術の開これらの条件を満たすため，筆者らはナノおよびメソ発ウェイトを上げ，3相界面やカーボン担体内部の超微小スケールのシミュレーションから得られた触媒層内の実細孔内の物理現象を正確に把握できるようにし，担体粒効的なプロトン伝導度と拡散係数を有する，プロトンと子の内部構造も含めた最適化に取り組み，性能と耐久の酸素の保存法則（1）式，（2）式を界面の酸素濃度を境界両面から見た定置用燃料電池の真の理想構造を明らかに条件として，あらかじめ解析的に解き，触媒層内での過していく所存である．

電圧の合計を出力とする過電圧関数を作成し，この過電参考文献

圧関数をマクロ化し，Fluentの燃料電池モジュールに組み込むとともに，MPL/GDLなどの各種多孔質材料の実効的な物質輸送パラメータを算出した．結果，第6図に示す
ように，筆者らのターゲットとなる全運転領域において，電流-電圧特性で実験との乖離が±10 mV以下になるこ76[1]

T. Mabuchi et al., “Effect of bound state of water on hydronium
ion mobility in hydrated Nafion using molecular dynamics
simulations,” J.Chem.Phys. vol. 141, no. 10, pp. 4904-4911,
2014.

77

先端研究特集：燃料電池マルチスケールシミュレーション[2]

T. Ohmori et al., “Anomaly of the temperature dependence of the

執筆者紹介

diffusion of oxygen in supercritical water,” J. Chem. Phys. vol.
116, no. 7, pp. 2680-2683, 2002.
[3]

P. Gode et al., “In-situ measurements of gas permeability in fuelcell membranes using a cylindrical microelectrode,” J.
Electroanal. Chem. vol. 518, no. 2, pp. 115-122, 2002.

[4]

A. Parthasarathy et al., “Investigations of the O2 Reduction
Reaction at the Platinum/Nafion® Interface Using a Solid‐StateElectrochemical Cell,” J. Electrochem. Soc. vol. 138, no. 4, pp.
916-921, 1991.

[5]

J. S. Chiou et al., “Gas permeation in a dry Nafion membrane,”Ind. Eng. Chem. Res. vol. 27, no. 11, pp. 2161-2164, 1988.

[6]

F.N. Büchi et al., “Microelectrode Investigation of Oxygen山本 恵一

Keiichi Yamamoto

先端研究本部
水素・エネルギープロジェクト室
Hydrogen And Energy Research Lab.,
Advanced Research Div.
藤田 悠介

Yusuke Fujita

先端研究本部
水素・エネルギープロジェクト室
Hydrogen And Energy Research Lab.,
Advanced Research Div.

Permeation in Perfluorinated Proton Exchange Membranes withDifferent Equivalent Weights,” J. Electrochem. Soc. vol. 143, no.
3, pp. 927-932, 1996.
[7]蔦原道久 , CFD最前線, 日本機械学会, 共立出版（株）, 東京,
2007．

[8]

安本栄一

Eiichi Yasumoto

アプライアンス社 技術本部
Corporate Engineering Div.,
Appliances Company

Li-Shi Luo et al., “Theory of the lattice Boltzmann method:Two-fluid model for binary mixtures,” Phys. Rev. E67, 036302,2003.

[9]E. P. Gross et al., “Model for Collision Processes in Gases:
Small-Amplitude Oscillations of Charged Two-Component
Systems,” Phys. Rev. vol.102 issue 3, pp. 593-604, 1956.

[10] Abhijit S Joshi et al., “Lattice Boltzmann method for continuum,

菅原 靖

Yasushi Sugawara

先端研究本部
水素・エネルギープロジェクト室
Hydrogen And Energy Research Lab.,
Advanced Research Div.

multi-component mass diffusion in complex 2D geometries,” J.
Phys. D: Appl. Phys. 40, pp. 2961-2971, 2007.

77