固体高分子形燃料電池シミュレータ P-Stack 次世代版

固体高分子形燃料電池シミュレータ P-Stack 次世代版の開発
米田雅一i
高山務ii
吉村英人iii
茂木春樹iv 高山糧v 仮屋夏樹vi
Development of Simulation Software for Polymer Electrolyte Fuel Cell
Masakazu YONEDA
Tsutomu TAKAYAMA
Hideto YOSHIMURA
Haruki MOTEGI Ryo TAKAYAMA Natsuki KARIYA
固体高分子形燃料電池は環境にクリーンなエネルギー革新的技術であり，次世代自動車の動力源やコジェ
ネレーションシステムとしての普及が期待されている．みずほ情報総研では数値シミュレーション技術を
活用した燃料電池の動作メカニズムの解析により運転条件・構造・部材特性に関して複雑なトレード・オ
フを持つ設計最適化の支援を目指し，固体高分子形燃料電池のシミュレータ P-Stack を開発・販売してい
る．ここでは，P-Stack を用いた単セルレベルでの検証，および 100 枚積層スタックの解析事例を紹介する．
(キーワード): 固体高分子形燃料電池（PEFC）
，燃料電池自動車（FCV）
，家庭用燃料電池，CAE，ソフトウェア，
セル・スタック，マクロモデル，気液二相流，膜電極複合体（MEA）
，ガス拡散層（GDL），設計解析
i
サイエンスソリューション部
デジタルエンジニアリングチーム
次長
ii
サイエンスソリューション部
デジタルエンジニアリングチーム
チーフコンサルタント
iii
サイエンスソリューション部
デジタルエンジニアリングチーム
コンサルタント
iv
サイエンスソリューション部
デジタルエンジニアリングチーム
コンサルタント
v
サイエンスソリューション部
デジタルエンジニアリングチーム
コンサルタント
vi
サイエンスソリューション部
デジタルエンジニアリングチーム
コンサルタント
6
博士（理学）
博士（工学）
はじめに
1
低減が必須の課題である．図 1 に燃料電池システム
の主な開発課題をまとめたものを示す．
燃料電池は，従来の内燃機関等に比べてエネルギ
燃料電池スタックのコスト低減のためには，一般
ー変換効率が高く，二酸化炭素排出の大幅な削減，
的に３つの方策が考えられる．具体的には，電極触
天然ガスや灯油，メタノール，バイオ燃料などの多
媒の貴金属使用量の低減や電解質膜のコスト低減等
様な燃料の使用が可能であり，石油代替の促進にも
による部材のコスト低減，燃料電池の高電流密度化
寄与する技術である．なかでも固体高分子形燃料電
による部材の使用量低減，そして，スタックの運転
池（PEFC：Polymer Electrolyte Fuel Cell）は出力電流
を制御する補機部品の簡素化，すなわちシステムの
密度が高く小型化・軽量化が可能，常温から発電で
簡素化である．
き起動停止性に優れる，作動温度が 80～90℃程度で
部材のコスト低減については，現在のスタックコ
耐熱温度の低い安価な材料が選択できる等の特徴か
ストで約 50%近くを占めると言われている白金触媒
ら，燃料電池自動車（FCV：Fuel Cell Vehicle）や家
の使用量を現状の 1/10 レベルまで低減することが
庭用コージェネレーションシステム（熱電併給）と
最も重要な課題の１つである．電解質膜はパーフル
しての実用化・本格普及が強く望まれているエネル
オロスルホン酸（PFSA）の膜の１つである Nafion 膜
ギー革新技術である．
が主流であるが，原材料が高く，かつ合成が煩雑で
2009 年に家庭用燃料電池「エネファーム」が市場
あることがコスト低減の壁となっている．近年では
投入され，2014 年度中には累計販売台数が 10 万台
炭化水素系電解質膜などの比較的コストが安価な電
を突破するなど順調に普及が進展しているなか，
解質材料が着目されているが，Nafion 膜と同等レベ
2014 年 12 月にはトヨタ自動車から一般向けで世界
ルのプロトン伝導性と耐久性向上が課題である．
初の FCV「MIRAI」の販売が開始された．本田技研
部材の使用量低減については，セルの発電面積や
工業も 2015 年度中の販売を表明するなど，水素エネ
積層した各セルへ燃料ガスを供給するマニホールド
ルギー利用社会の新しい幕開けが到来しつつある．
断面積を小さくして，スタック全体をコンパクト化
当社では，数値モデリングとシミュレーションに
することが課題である．しかし，局所的な発電集中
よる燃料電池開発の支援を目的として，燃料電池の
や水分不足あるいは水分過多は発電特性と部材の耐
発電特性と内部現象を予測するソフトウェア（固体
久性に影響を及ぼすため，セル内部の発電分布や水
®
高分子形燃料電池シミュレータ P-Stack ）を開発し，
分分布，セル間の発電性能の均一化を維持するため
ソフトウェアの年間使用許諾ライセンス提供および
の設計が重要となる．スタックの設計には，燃料ガ
P-Stack を用いた燃料電池セルの性能解析や燃料電
スを供給する流路構造，発電によるスタック温度上
池システムに係わる物理モデル開発，受託解析コン
昇を抑制するための冷却材を供給する流路構造など，
サルティングを行っている．
広範な形状設計パラメータの検討が必要となるため，
P-Stack はこれまで実機セルの構造・運転条件違い
CAE（Computer Aided Engineering）の設計解析への
に対する発電特性および反応と熱物質輸送の連成効
活用が有効である．
果を解析することが一つの目標であったが，数値解
最後にシステムの簡素化については，燃料電池ス
析モデルの開発強化により，最終的には 400～500
タックが安定に作動する各種条件（出力密度，温度，
枚セルを積層した FCV の燃料電池スタックまで発
圧力，燃料ガス流量等）の範囲を拡大し，補機部品
電特性解析の対象を拡張して開発を進めている．本
の汎用品への切り替えや部品点数を削減してコスト
稿では P-Stack の導入メリットと特徴，数値解析モ
低減を図ることが課題である．一例として，PEFC
デルの概要と検証結果，実機セルおよびスタックレ
では 80～90℃程度の温度を上限として燃料電池ス
ベルの解析事例について紹介する．
タックを作動させるが，外気温度との差が大きいた
めに放熱条件が従来の内燃機関を用いた自動車と比
燃料電池開発における技術的課題
2
べて厳しく，ラジエータシステムが大型になる．こ
のため，更に高温で燃料電池スタックを作動するこ
2.1
コスト低減に向けた解析技術の有効性
とが要求されるが，水の沸点（100℃）を超える条件
固体高分子形燃料電池の普及拡大に向けて，燃料
では発電による生成水の持ち去りが急激に大きくな
電池の発電性能と耐久性能の向上の両立と，コスト
るため，電解質内の水分を保持することが非常に難
7
しくなる．さらに，水分制御の技術が向上すれば，
影響因子
燃料ガスを加湿せずに発電による生成水のみで電解


質膜のプロトン伝導性を維持することも可能となり，
発電性能
評価項目
加湿器の除去等，大幅なシステムの簡素化を図るこ
とができる．FCV の実運転モードにおけるスタック


の発電特性と内部の反応・水分・ガス濃度の変動を
予測することが可能となればシステム設計にも有効

である．

発電性能向上
高電流密度化
耐久性能
評価項目
コンパクト化


材料コスト低減
貴金属使用量低減
構造条件：MEA諸元、流路パターン、流路構造（断面形状、
表面濡れ性）、マニホールド構造、等
作動条件：ガス・冷却水流量、圧力、加湿、温度、電流密
度、負荷変動パターン、等
スタックの発電分布・含水分布の均一性、セル間流量バラ
ンスの設計基準を満たすための流路形状・マニホールド
形状の評価
低ストイキ（水素デッドエンドを含む）、高温（100℃以上）
作動時のスタック内物質移動（水分・ガス濃度）バランスと
発電分布の均一性の評価
Pt量削減時のMEA特性の変化に応じたセルおよびスタッ
ク内発電分布への影響の評価
低温起動時のスタック昇温時間を満たすための発熱・伝
熱の評価、等
流路・マニホールド形状に対する起動・停止時の水素置換
状況の評価およびカーボン劣化量の評価
負荷変動時に発生する乾湿変動、過電圧変動、温度変動
が大きい領域の推定、等
低ストイキ化
トレード・オフ両立
水素循環レス
燃料電池シミュレータ P-Stack の概要
3
システム簡素化
高温・加湿器レス
低温起動性向上
耐久性能向上
3.1
燃料電池の構造と内部現象
燃料電池シミュレータは電池内部を支配する熱・
耐久性能向上
物質の輸送と電気化学反応の連成効果を数値モデル
触媒（起動・停止、負荷変動）
化することによって，燃料電池の構造や作動条件に
電解質（高電位保持、乾湿変動）
おける電流密度分布，ガス濃度分布，水分分布等の
セパレータ・シール材強度、熱耐性
内部状態と分極特性（I-V 特性）等の発電特性を予
図 1 燃料電池システムの開発課題
測するものであるが，考慮すべき動作メカニズムは
広範である．図 2 に PEFC の構造とそれを構成する
2.2 燃料電池シミュレータの導入メリット
前節で説明したように，燃料電池のコスト低減に
向けた開発・設計にシミュレーションを活用するこ
とが有効である．具体的には，
部材，各部材における現象を示す．PEFC の発電性
能に影響を与える現象は，発電部の心臓である膜電
極複合体（MEA : Membrane Electrode Assembly）の
電解質膜，触媒層に関するミクロ・ナノレベルの現



象から，ガス拡散層（GDL : Gas Diffusion Layer）の
実験では「見えない」内部現象理解による設計
指針の策定支援
新規開発のセル・スタック性能の事前予測によ
る試作コスト，開発工程の短縮
実験では部材劣化の懸念がある限界条件での
性能予測
カーボン繊維のメゾレベルの現象，さらにセパレー
タの流路パターンと伝熱・流動，燃料電池スタック
全体のガス・冷却水の供給バランス，発電による発
熱と放熱バランス等のマクロレベルの現象まで，幅
広いスケールの現象を含んでいる．同時に，触媒層
における電気化学反応，電解質内の水分とプロトン
を可能とし，表 1 に示すアウトプットに対して高速
の輸送，GDL および流路内のガスと生成水の二相流
かつ設計に要求される精度で解析できることを目標
動を含む物質輸送と熱輸送など，内部では複雑な現
に開発を進めている．
象が連成している．このため，PEFC の挙動をシミ
ュレーションするためには，ミクロ・ナノレベルか
表 1 P-Stack を用いた設計・評価項目事例
らマクロレベルまでの各スケールに適切なモデリン
グを行い，それらを連携させるシミュレーション技
術の枠組みが必要である．
8
図 3 100 枚積層スタック温度分布
PEFC 内の熱伝達および物質輸送と電気化学反応
を連成して計算する方法としては，通常の 3 次元流
体解析（CFD : Computational Fluid Dynamics）に基づ
いた数値モデル 1-5)が考えられるが，実機セルの発電
特性をシミュレーションするためには膨大な計算資
源と計算時間を必要とし，さらに前述のように MEA
の基本特性（作動条件に対する I-V 性能）をシミュ
(a) 燃料電池スタック構造
レーションによって導出するためには，ナノスケー
V
e-
燃料極（アノード）：H2 → 2H+ + 2e -
空気極（カソード）：O 2 + 4H+ + 4e -→2H2 O
H2
e-
ルレベルの構造と輸送現象を理解し，部材要素ごと
の物性値を連続体モデルとして構築する必要がある．
O2
H+
このため，P-Stack では，発電の要素特性を担う
H2
H2O O
2
MEA 単体の基本特性については小型セルの計測デ
Drag
H2
ｱﾉｰﾄﾞ
ｾﾊﾟﾚｰﾀ・流路
ｱﾉｰﾄﾞ
GDL
ｱﾉｰﾄﾞ
MPL
ータを活用し，セルの性能を決定づける MEA の発
O2
Diffusion
ｱﾉｰﾄﾞ
電解質膜ｶｿｰﾄﾞ
触媒層
触媒層
ｶｿｰﾄﾞ
MPL
ｶｿｰﾄﾞ
GDL
電特性や GDL 内の物質輸送特性に関するマクロモ
ｶｿｰﾄﾞ
ｾﾊﾟﾚｰﾀ・流路
デル（集中定数系モデル，3.3 節を参照）のパラメー
MEA
セル
タを同定し，流路・マニホールド内の流量と圧力損
失を決定づける流動パラメータは実形状に対して
(b) セルを構成する部材と輸送現象
部材名
ｾﾊﾟﾚｰﾀ・流路
GDL
MPL
触媒層
電解質膜
寸法
～1.0mm
～300μm
～100μm
～10μm
～50μm
機能
・ｶﾞｽ供給（流路）
・生成水排出（流路）
・ｾﾙ間の熱・電気ｺﾈ
ｸﾀｰ
・反応ｶﾞｽ拡散
・生成水排出
・熱、電気伝導
・反応ｶﾞｽ拡散
・生成水排出
・熱、電気伝導
・電極表面反応
・反応ｶﾞｽ拡散
・生成水排出（空孔）
・熱、電気伝導（電極）
・ﾌﾟﾛﾄﾝ伝導（ｱｲｵﾉﾏｰ）
・ﾌﾟﾛﾄﾝ伝導
・水分輸送
（濃度拡散、電気浸透）
・ｸﾛｽｵｰﾊﾞｰ抑止
・ｶｰﾎﾞﾝ含浸樹脂＋
表面撥水処理
・SUS系ﾌﾟﾚｽ形状
・ｶｰﾎﾞﾝ繊維
（繊維径10μm）
・ﾍﾟｰﾊﾟｰ/ｸﾛｽ
・PTFE処理
・ｶｰﾎﾞﾝﾌﾞﾗｯｸ
（粒子径0.1μm）
・PTFE処理
・Pt担持ｶｰﾎﾞﾝﾌﾞﾗｯｸ
（粒子径0.1μm）
・ｱｲｵﾉﾏｰ（PFSA系、等）
・Nafion（PFSA系）
・炭化水素系
構造
材質
CFD を実施し，内部の速度・圧力分布の結果から同
定する．
さらに，触媒層の微細構造等に起因した高負荷領
域での限界電流の発生，アイオノマーのプロトン輸
送抵抗や細孔内の酸素輸送抵抗による濃度過電圧の
発生等をモデル化して MEA の基本的な挙動を再現
(c) 各部材の機能と構造・材質
できるように拡張 6)している．
図 4 に P-Stack における計算格子のイメージを
図 2 PEFC の構造（部材構成）と内部現象
CFD の計算格子と比較して示す．P-Stack では触媒
3.2
層（CL：Catalyst Layer）を表面モデル（0 次元モデ
機能・特徴
P-Stack は，発電面積が数百 cm のセルおよびスタ
ル）
，流路についてはセル面方向の 2 次元配管モデル
ック（FCV 動力用 PEFC では数百枚積層，家庭用
で扱う．一方，GDL は 3 次元的に扱い，リブ直下の
PEFC では数十枚積層）の作動条件と構造条件に対
酸素濃度低下と含水量の増加，流路間圧力差に起因
する電流－電圧特性（I-V 特性）と負荷変動特性，
する GDL 面方向のガスのパスカット等の解析にも
内部の電流密度分布，ガス濃度・水分・温度分布等
対応している．
2
を解析するシミュレータである．図 3 に出力電流密
度 0.2A/cm2 における 100 枚セル積層スタックの温度
流路
分布を示す．
水素
空気
GDL+MPL
353 K
CL
PEM
冷却水
348
(a) CFD に基づく計算格子
流路
GDL+MPL
(CL)
PEM
9
おける含水量との関係を計測データから決定し，
（b）P-Stack における計算格子


以上のプロセスにより，MEA 諸元や流路パターン
(5)
抵抗過電圧は，電解質厚さ方向で積分する．
等の部材・構造条件，作動条件の違いによる実機セ
 r ( m)  
ルあるいはスタックの発電特性と内部分布の連成効
i
t( m )
0
果を高速に解析することが可能となる．
 ( m) ( z )
(6)
dz
電気浸透および濃度拡散による極間の水移動量は以
下の式を適用し，電解質膜と GDL との界面は水分
数値解析モデル
3.3
 1
1 


 303 T( m) 


 (effm)   w( m) exp1268
図 4 燃料電池モデル体系と格子解像度のイメージ
輸送抵抗 k weff を用いて境界条件としてモデル化する．
P-Stack における熱・物質の輸送と電気化学反応の


 d
H 2O
 d z M ( m)  0


dry
 M H 2O  neff i  ( m ) D eff d w( m )
d
w( m )
 (m)
F EW
dz

連成効果を計算するための MEA の起電力モデル，
電解質膜内の物質輸送モデル，流路・マニホールド
および GDL 内の熱・物質輸送モデル 7)について以下
に示す．

p
H 2O 
kweff( a )  s a w( m)  cGDL
(a) 


 RT


H 2O 
 i  k eff  ps a 
w( c ) 
w( m )  cGDL ( c ) 
 2 F
RT



3.3.1 起電力モデル
M *H 2O
セル電圧と電流密度の関係は，開回路電圧，両極
の電気化学反応による活性化過電圧，電解質膜中の
抵抗過電圧，および触媒層の微細構造等に起因した
濃度過電圧（カソード）である．
Vcell  Voc  a( a)  a(c)  r ( m)  c(c)
3.3.2
(7)



(8)
GDL 内熱・物質輸送モデル
(1)
流路および GDL 内のガスと液水輸送の基本モデ
起電力の基本モデルは以下の Butler-Volmer 式であ
ルは，二相流の圧力損失式と質量保存式である．
GDL はメゾスケールのカーボン繊維構造に起因す
る．
るガスと液水の輸送挙動について，毛管圧力，相対
   a( s ) F

  a( s ) F

i  i0( s ) exp
 a ( s )   exp 
 a ( s ) 

 RT( m)

  RT( m)



浸透率および有効ガス拡散係数をパラメータとした
多孔質二相理論モデルで表現している．
圧力損失式（Darcy 則）：
(2)

i0( s )  i0ref( s ) (effs) 


p(Pts)
p(ref
s)




 (s)
 G( s )
exp
 R
 1
1


 T( m) T ref


K eff
1  s 3 p g
 ug  
g


K eff 3

u


s p g  pc ( s )
 l
l






(3)


i
p(Pts)  p( s ) 1  eff

 ilim ( p( s ) ) 
(9)

質量保存式：
(4)

    g 1  s      g 1  s u g
 t

H O
k k
  w m p  w 2 m v
 k
 
H O
   l s      l s ul    w 2 m v
  t

式(1)，式(2)で表される Butler-Volmer 式と交換電流
密度は素反応レベルの物理化学的な性格を持つが，
実際の触媒層はカーボン粒子の凝集・分散構造，Pt
の粒径分布，アイオノマー被覆構造等，ナノレベル
で複雑な構造を持っている．ここではこれらの複雑
な構造に起因する触媒特性をマクロモデルで扱うた



(10)
ここで，m kp は電気化学反応による成分 k の生成・消
滅量，m v は蒸発・凝縮による相間質量移動量である．
めに低加湿状態での触媒利用率の低下，式(4)で表さ
れる反応ガス低下による限界電流密度の低下をパラ
 i 
m kp     k 
n F 
メータとして拡張し，これを計測データにより決定
する．また電解質内のプロトン伝導度は基準温度に
10
(11)
p 

 v  kveff ( s) wH 2O  c H 2O  s 
m
RT 

scr 
(12)
 d cr3
(18)
6vo
濃度保存式：
d cr については，可視化セル試験による文献 8)も見ら
 k
れるが，ここでは VOF 法や粒子法 9, 10)等の直接二相

c 1  s     c k 1  s u g  Dgk ,eff c k
t
流解法を用いた流路と GDL 界面の液水挙動に関す
(13)
k
る二相流解析の結果 11)からモデル化した．
 m p  m v
図 5 に，カソード流路内（プレス成形セパレータ
ここで Dgk , eff は多孔質内の有効ガス拡散係数であり，流路）のガス流速に対する GDL 界面の限界体積率
との相関，液水体積率と圧損増倍との相関をそれぞ
GDL の空孔度  と屈曲度  ，および液水滞留による
フラッディングを仮定した以下のモデルを適用した．れ示す．流速が増加すると疎水性の流路表面である
ほど限界径が小さくなって排水性が向上する．一方，

Dgk ,eff  Dgk   1  s 
(14)
流速が遅い領域では親水性の流路表面のほうが
エネルギー保存式：
GDL 表面に排出された水滴を吸い上げる影響によ


1    

t

 s hs  



 g 1  s hg  l shl
t



って，排水性は高くなっている．圧力損失について
は疎水性の流路表面では水滴が流路を閉塞して圧損

    g 1  s hg ug  l s hl ul     T  q( m)
増倍が大きくなるのに対し，親水性の流路表面では
(15)
水滴が液膜状となるために圧損増倍は小さくなる．
eff
0.40
Maximum Saturation / -
また，図 3 に示すように，流路－GDL 界面における
熱・物質輸送フラックスは以下のようにモデル化さ
れる．なお，Sherwood 数および Nusselt 数は流路内
が層流流れを前提とした値を用いる．
 k
Dgk ,eff k
k
cchannel  cGDL
 M *  Sh

h

eff k

k
q

Nu
Tchannel  TGDL
 *
h





(16)
Hydrophobic, θ=120
Hydrophilic, θ=60
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.00
0
5
10
15
20
Gas Velocity / m/s
(a)
GDL 界面の液水限界率と流速との相関
3.3.3 流路内熱・物質輸送モデル
Two-Phae Pressure Increase / -
流路内の熱・物質輸送モデルは式(9)～(15)に対し
て空孔度を 1.0 として扱う．だたし，圧力損失式は
水力学的等価直径 d e を用いて流路の層流粘性抵抗
をモデル化した以下の式を適用した．

d e2
p g
 ug  
32 g TP ( s )


0
: s  scr




2
 ul   d e sp : s  s
g
cr


 32l

(17)
3.0
Hydrophobic, θ=120
Hydrophilic, θ=60
2.5
2.0
1.5
1.0
0.5
0.0
0.00
0.05
0.10
0.15
0.20
0.25
Liquid Water Saturation /-
(b) 液水率と圧損増倍係数との相関
ここで， TP は液水存在時における圧損増倍係数で
図 5 流路二相流モデルパラメータ例
あり，また，液水がガス流れによって移動を開始す
る限界体積率 scr は，発電時に GDL から排出された
3.4
液水の限界成長径 d cr （detachment diameter）から決
P-Stack を用いた解析プロセス
P-Stack を用いた解析プロセス（手順）を図 6 に示
定した．
す．最初に，セルユニットを構成する 3D-CAD デー
11
タ（セパレータ，GDL，MEA およびシール材）を
インポートし，図 7 に示すように，P-Stack メッシュ
データおよび熱流動パラメータ同定用の CFD メッ
シュデータを自動的に生成する．
P-Stack のメッシュ解像度はセパレータの形状に
もよるが，最小で 0.1mm 程度（エンボス・整流板周
(a) P-Stack メッシュ（単セル全体）
りなど）
，最大で 5mm 程度（流路流れ方向の間隔）
として計算負荷を軽減している．一方，CFD のメッ
シュ解像度は P-Stack の約 1/10 であり，セル１枚あ
たりのメッシュ規模は最大で P-Stack の 1,000 倍とな
る．
次に，セルユニット単位で CFD を実行して内部の
(b) CFD メッシュ
流速分布と圧力分布を計算し，流路領域における式
（流路－バッファ接続・マニホールド領域の拡大）
(17)の水力学的等価直径を算出する．この段階で，
セル・スタック構造に対する流動評価（流路間およ
図 7 メッシュ解像度イメージ
びセル間流量分布および圧力損失）が可能となる．
さらに，実機セル・スタックの発電性能を解析す
るために，実機で想定される酸素濃度，加湿度およ
モデル検証と解析事例の紹介
4
び温度条件の作動範囲に対する小型セルの電流－電
圧（I-V）特性と電流－抵抗（I-R）特性を計測し，
ここでは，セルレベルでの数値計算モデルの基本
その結果を再現するように MEA の電気化学反応パ
検証結果について説明し，実測データを再現するモ
ラメータや GDL 内の物質輸送特性に関するモデル
デルパラメータを用いて，実機サイズセルおよび
パラメータを同定する．
100 枚積層スタックに対する発電特性の解析を実施
した事例について紹介する．
これらのモデルパラメータを用いて，実機セルお
よび実機スタックの発電性能の解析を行うことがで
きる．
4.1
数値計算モデルの基本検証
ここでは、単セルの実測データを用いて，前章で
説明した数値解析モデルの基本検証を実施した．ま
CADデータ
ず，入口酸素濃度（21, 10, 5%）と入口加湿度違い（100,
メッシュ生成（自動処理）
70, 40%RH）の 9 条件に対する発電面積 25cm2 の単
メッシュデータ（P-Stack＋CFD）
セル（JARI）の I-V 特性を計測し，その結果と定量
的に一致するように式(1)～(18)の反応・物質輸送パ
CFD計算（単セル）
ラメータを適合した 7)．
熱流動パラメータ同定
なお，セル面内の発電分布の均一性の確保と流路
セル・スタック熱流動解析
実験データ
内の液水滞留を抑制するために，水素と空気の流量
圧損、流配、熱バランス出力
は利用率が 10%以下となるように条件設定した．ま
MEAパラメータ適合
たセル温度についてはほぼ 80℃固定となるように
発電特性パラメータ同定
調整した．
セル・スタック発電性能解析
表 2 に MEA（触媒層・電解質膜）および GDL の
モデルパラメータの適合結果を，図 8 に各条件にお
発電特性・内部分布出力
ける I-V 特性の比較をそれぞれ示す．全ての条件に
図 6 P-Stack による解析プロセス
おいてほぼ定量的に一致しており，これらの結果か
ら P-Stack モデルが実測 I-V データの再現性を有して
いることが判断できる．
12
表 2 反応・輸送モデルパラメータの適合結果
パラメータ
1.0
適合値
RH40_O2_21%(cal.)
measured
RH40_O2_10%(cal.)
measured
RH40_O2_5%(cal.)
measured
0.9
触媒層
界面水分輸送抵抗
アノード反応
カソード反応
i0( a )  1.0  10 4
i0( c )  3.2  10 1
 (a )
 (c)
 2.0
0.8
Voltage / V
電気化学反応
パラメータ
（Butler-Volmer式）
 0.81
 ( a)  0.5
 (c)  0.48
G( a )  0.0
G(c)  4.5 10 4
p Href2  1.01325  105
pOref2  1.01325  105
k weff( a )  2.0  10 5
k weff( c )  4.0  10 5
0.6
0.5
0.4
0.3
 1.5 pO2
eff
10
ilim(
a )  1.0  10
eff
ilim(
c)
0.7
0.2
0
0.0
電解質膜
電気浸透係数
nd  1.0
Dweff( m )


6000
0.6
8000
0.8
10000
1.0
12000
1.2
(c) 低加湿条件（RH=40%）における酸素依存性
 3.1  10 7  e 0.28  1  e 2346 / T : 0    3

 4.17  10 8  1  161e   e  2346 / T :   3

4000
0.4
Current Density / A cm-2
 0.514  0.326 e
プロトン伝導度
水分拡散係数
2000
0.2
1268 1 / 3031 / T 
 (effm )

図 8 電流－電圧特性の実測値と計算値との比較
GDL
有効ガス拡散係数
Dgeff  Dg  1  s 3.0
透過係数
K eff  8.0 10 12
次に，カソード流路の表面に酸素濃度検出試薬
Platinum-tetrakis Pentafluorophenyl Porphyrine (PtPP)
を塗布した可視化セル（I-V 特性と同じ JARI セルベ
ース）を用いて，流量（300mL/min）固定時の出力
電流密度の増加に対する流路内酸素分布を計測 12)
し，計算結果と比較した．なお，入口加湿度 60%RH，
セル温度 80℃で一定とした．
図 9 に各出力電流密度（0.22, 0.43, 0.58A/cm2）に
おける発電分布の変化を，図 10 に各出力電流密度に
おけるカソード流路内酸素分圧分布の変化をそれぞ
れ示す．出力電流密度が低い条件では，発電分布が
均一となるために流路内の酸素分圧は入口からほぼ
一定の傾きで減少する．一方，出力電流密度が増加
して酸素利用率が高くなると，発電分布がカソード
入口側に集中してその領域での酸素消費量が増加す
るために，流路内の酸素分圧は上流域での減少度が
大きくなり，分布は下に凸の様相となる．
P-Stack による計算結果は実測データを定量的に
再現できており，出力電流密度違いによる反応と物
質輸送の連成効果を表現する数値計算モデルとして
妥当であるといえる．
1.0
RH100_O2_21%(cal.)
measured
RH100_O2_10%(cal.)
measured
RH100_O2_5%(cal.)
measured
0.9
Voltage / V
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0
0.0
2000
0.2
4000
0.4
6000
0.6
8000
0.8
10000
1.0
12000
1.2
Current Density / A cm-2
(a) 高加湿条件（RH=100%）における酸素依存性
1.0
RH70_O2_21%(cal.)
measured
RH70_O2_10%(cal.)
measured
RH70_O2_5%(cal.)
measured
0.9
Voltage / V
0.8
0.7
0.215
0.222
0.365
0.449
0.183
0.766 A/cm2
0.6
0.5
0.4
0.3
0.2
0
0.0
2000
0.2
4000
0.4
6000
0.6
8000
0.8
10000
1.0
UO2=30% (0.22 A/cm2)
12000
1.2
Current Density / A cm-2
60% (0.43 A/cm2)
80% (0.58 A/cm2)
図 9 出力電流密度に対する内部発電分布の変化
(b) 中加湿条件（RH=70%）における酸素依存性
13
0.0
タについては，セパレータが SUS 系材料のプレス形
15.0 kPa
状のタイプであるため，式(17)の親水性表面におけ
る相関式を適用した．
(a) アノード流路（燃料極）
UO2=30% (0.22 A/cm2)
60% (0.43 A/cm2)
80% (0.58 A/cm2)
(a) 酸素分圧分布のコンター
（上段：可視化セル結果／下段：計算結果）
(b) カソード流路（空気極）
(c) 冷却水流路
図 11 実機セルの流路パターン
表 3 セパレータの流路構造の違い
項目
設定値
セルサイズ
 全体サイズ：310mm×150mm
（31本並行流路）
 MEA面積: 260cm2
セパレータ
 材質：プレス成形セパレータ（SUS）
 流路深さ：
(a) An=0.5mm, Ca=0.5mm, Coolant=0.5mm
(b) An=0.4mm, Ca=0.4mm, Coolant=0.4mm
 流路ピッチ：3.9mm
 リブ幅：1.5mm
 バッファ（分岐・合流）領域深さ：
(a) An=0.3mm, Ca=0.3mm, Coolant=0.2mm
(b) An=0.2mm, Ca=0.2mm, Coolant=0.2mm
(b) カソード入口からの距離に対する酸素分圧分布
図 10 出力電流密度に対する酸素分圧分布の変化
4.2
4.2.1
実機サイズセルの解析事例
解析条件
ここでは，実機サイズセルの構造条件違いに対す
る発電特性および内部分布への影響を解析した事例
を紹介する．
図 11 に各流路のパターン（並行流路）を示す．発
電面積（点線領域）は 260cm2 とし，表 3 に示すよう
4.2.2 解析結果
に流路の薄厚化（セルユニットのコンパクト化，流
流路構造違いによる負荷変動時のセル電圧変化を
路深さ：セル"a"＞セル"b"）による発電性能への影
図 12 に示す．また，出力電流密度 0.25A/cm2 および
響を解析した．
運転条件は入口利用率
（アノード 70%， 1.0A/cm2 における発電分布の変化を図 13 に，MEA
カソード 50%）
，入口加湿度 RH=50%とし，出力電
2
面内の 6 点（カソード入口側 3 点，カソード出口側
2
流密度 0.1A/cm から 0.25, 0.5, 1.0A/cm と負荷変動
3 点）の直下における電解質膜内の水蒸気活量の時
させたときのセルの発電特性および内部状態の非定
間変化を図 14 にそれぞれ示す．
常挙動を解析した．
出力電流密度を増加させたときの電圧変化は，セ
なお，MEA と GDL のモデルパラメータは表 2 に
ル"b"のほうが低下する．これはセル"b"ではセル"a"
示した適合結果を用い，流路二相流モデルパラメー
と比べると，流路深さ以上にバッファ領域深さが大
14
1.00
▲33%）
，流路よりバッファ領域の圧力損失が大きく
0.95
Water Actitivy in PEM / -
きく減少するため（流路深さ▲20%，バッファ深さ
なり，各流路への流量分配が不均一になるためであ
る．
具体的には，セル両端の流路（図 12 の A・C）を
通過する流量が増加し，セル中央の流路（図 12 の B）
を通過する流量が減少するため，セル両端の入口領
域では水蒸気活量が低下して電解質膜抵抗が増加，
一方，セル中央の出口領域では水蒸気活量が過多と
In-A
In-B
In-C
In-A
In-B
In-C
0.90
0.85
Cell-(b)
Cell-(a)
0.80
0.75
0.70
0.65
0.60
0.55
0.50
なってフラッディングによる酸素拡散の阻害が発生
0
5
10
15
して発電分布が面内で不均一になる．
Water Activity in PEM / -
た，電解質膜内の含水状態の乾湿サイクルをシミュ
レーションすることにより，負荷変動による電解質
膜の機械的劣化が発生しやすいポイントの予測も可
1.0
0.75
0.8
0.70
Current Density
0.6
0.65
Voltage : Cell (a)
Voltage : Cell (b)
0.4
0.60
0.2
0.55
0.0
5
10
15
20
25
30
1.15
1.10
1.05
1.00
Out-A
Out-B
Out-C
Out-A
Out-B
Out-C
0.95
0.90
Cell-(b)
Cell-(a)
0.80
0
5
10
15
20
25
30
35
Time / sec
(b) カソード出口領域
図 14 電解質膜内水蒸気活量の時間変化
4.3
0.50
0
1.20
0.85
Cell Voltage / V
Current Density / A cm-2
能となる．
0.80
35
1.25
一性は増加して電圧が低下することがわかった．ま
0.85
30
1.30
流密度増大に伴い，面内の局所電流密度分布の不均
1.2
25
(a) カソード入口領域
以上から，セルを薄厚化することによって出力電
1.4
20
Time / sec
35
セル 100 枚積層スタックの解析事例
4.3.1 解析条件
Time / sec
図 11 に示した流路パターンで構成されるセルを 1
図 12 負荷変動時のセル電圧変化
-30
セルのユニットとして 100 枚積層したスタックにつ
いて，各セルへ水素，空気および冷却水を供給する
+30 %
マニホールド断面積の違い，および出力電流密度の
基準セル”a”
In-C
Out-C
違いによる発電特性への影響を解析した事例を紹介
In-B
Out-B
する．マニホールド断面については基準幅を 20mm
In-A
Out-A
として 15mm と 10mm の 3 形状を対象とし，出力電
流密度については低負荷の 0.2A/cm2 の場合と高負荷
薄厚セル”b”
の 1.0A/cm2 の場合の作動条件を対象とした．
In-C
Out-C
In-B
Out-B
In-A
Out-A
出力電流密度
0.25A/cm2
図 15 にセル積層構造とガス・冷却水の流れ方向を
示す．今回、
解析対象としたスタックは 1 セル
（MEA）
に対して 1 冷却の積層構造であり，各セルにおいて
出力電流密度
アノードとカソードは対向流，冷却水はカソードと
1.0A/cm2
同じ方向の流れと設定した．また，冷却水マニホー
図 13 発電分布の変化
ルド入口側のセル番号を No.1 として積層方向の番
号付けを行った．
15
-30
空気
冷却水
+30 %
水素
カソード流路
MEA＋GDL
アノード流路
冷却水流路
図 15 セル積層構造
0.780
8.0
0.775
6.0
0.770
4.0
0.765
2.0
0.760
0.0
Manifold Width=20mm
Manifold Width=15mm
Manifold Width=10mm
Manifold Width=20mm
Manifold Width=15mm
Manifold Width=10mm
0.755
0.750
0.745
-2.0
-4.0
0.740
0
10
20
30
40
50
60
70
80
90
図 17 エンドセルにおける発電分布
（上段：マニホールド幅 20mm／下段：10mm）
4.0
27.0 kPa
図 18 エンドセルにおけるカソード酸素分圧分布
（上段：マニホールド幅 20mm／下段：10mm）
Cathode Flow Rate / %
Cell Voltage / V
4.3.2 解析結果
図 16 にマニホールド断面積の違いに対する出力
電流密度 0.2A/cm2 におけるセル積層方向の電圧分布
とカソード流量分布を示す．また，図 17 と図 18 に
エンドセル（セル No.1：空気の流入方向に対する末
端セル）における発電分布とカソード酸素分圧分布
を示す．
マニホールドが 20mm 幅の場合は各セルの流量が
±2%以下で分布しているが，マニホールド幅が小さ
くなるほど，マニホールド内の圧力損失が増加する
ために各セルへの供給される空気流量の分布が不均
一となり，10mm 幅では 20mm 幅に比べて分布が 3
倍程度に拡大する．
特にエンドセル付近に流入する空気流量が低下し
ており，同エンドセル内部ではセル中央の下流域で
酸素不足が顕著に現れ，発電分布が不均一となって
いる．また，エンドセルではスタック端面からの放
熱が大きいために影響で液水滞留量が増加し，セル
内の酸素濃度低下を拡大させていると考えられる．
次に，マニホールド断面を基準幅（20mm）とし
て出力電流密度 0.2A/cm2 と 1.0A/cm2 に対する発電
特性と内部分布への影響を解析した．
図 19 に出力電流密度の違いによる各セルの電圧
-6.0
分布を，図 20 に出力電流密度 1.0A/cm2 における電
-8.0
解質膜内の抵抗過電圧，およびカソード触媒層の反
100
応過電圧分布をそれぞれ示す．なお，各セルの過電
Cell Number / -
圧は MEA の平均値である．また，図 21 と図 22 に
図 16 マニホールド断面積に対する電圧・流量分布の比較
（出力電流密度 0.2A/cm2）
出力電流密度 1.0A/cm2 におけるセンターセル（セル
No.50）とエンドセル（セル No.1）における MEA（シ
ール材を含む）層の温度分布およびカソード GDL
16
内に滞留する液相体積率分布をそれぞれ示す．
348
362 K
2
図 19 に見られるように，
出力電流密度が 0.2A/cm
の条件ではセル間の電圧差にバラツキが無いが，出
力電流密度が高くなるとエンドセルの発電性能が相
対的に低下する．さらにはエンドセル付近からセン
ターセル前後まではセル間の電圧差が小さくなって
いる．これは図 21 に見られるように，出力電流密度
が高くなると発電による発熱量が増加してセンター
セルの温度が上昇し，電解質膜抵抗の増大に繋がっ
ていることが要因である．一方，エンドセルでは液
水滞留量の増加（図 22）によって酸素拡散抵抗が増
加し，反応過電圧の増大を引き起こしている（図 20）．
このように高負荷条件においてはセル間の温度差に
起因して，セル間の発電性能差およびスタック全体
の発電特性に大きな影響を及ぼす可能性があること
図 21 出力電流密度 1.0A/cm2 における
MEA・シール材の温度分布
（上段：センターセル／下段：エンドセル）
が示唆された．
0.780
0.520
0.775
0.0
0.515
1.0 A/cm2
0.770
0.510
0.765
0.505
0.760
0.500
0.755
0.495
0.750
0.490
0.745
0.485
0.740
20.0 %
Cell Voltage @1.0A cm-2 / V
Cell Voltage @0.2A cm-2 / V
0.2 A/cm2
0.480
0
10
20
30
40
50
60
70
80
90
100
Cell Number / -
図 19 出力電流密度に対する電圧分布の比較
0.720
0.05
ORR in the Cathode CL / V
Ohmic Loss in the PEM
0.710
0.04
0.705
0.700
0.03
0.695
0.690
0.02
0.685
0.680
0.01
Ohmic Loss in the PEM / V
ORR in the CCL
0.715
図 22 出力電流密度 1.0A/cm2 における
カソード GDL 内液水分布
（上段：センターセル／下段：エンドセル）
5
0.675
おわりに
固体高分子形燃料電池シミュレータ P-Stack の次
0.670
0.00
0
10
20
30
40
50
60
70
80
90
世代版として燃料電池スタックの開発・設計に向け
100
Cell Number / -
た導入メリットと特徴，数値解析モデルの概要と検
図 20 出力電流密度 1.0A/cm2 における
膜抵抗過電圧とカソード活性化過電圧分布の比較
証結果，実機セルおよびスタックレベルの解析事例
について紹介した．
また，P-Stack では今回紹介した解析事例以外にも，
起動時の水素置換，発電時の局所スタベーションな
ど，電極の耐久性能に関する解析も可能である．
17
今後，3.4 節で説明した 3D-CAD インポートおよ
びメッシュ生成機能を含む次世代版のリリースを
2015 年 8 月末に予定しているが，それ以後について
も，燃料電池の開発・設計現場で有効となる 3D-CAE
との連携機能の強化（流路内二相流の高度化，構造
解析による部材変形と流路断面・接触抵抗の変化，
MEA 乾湿変動データによる電解質膜の機械的耐久
性評価，流路内二相流の高度化，など）を図る計画
である．
使用記号
i
i0
a
：電流密度 [A/m2]
：交換電流密度 [A/m2]
：活性化過電圧 [V]
ps
：飽和蒸気圧（温度関数） [Pa]
D
：ガス拡散係数 [m2/s]
h
：エンタルピー [J/kg]

q
d
 TP
vo
Sh
Nu
：熱伝導率 [J/m/K]
ref
Pt
eff
(s )
(m)
(g )
(l )
k
：基準状態参照値
：反応による発熱量 [J/m3/s]
：水力学的等価直径 [m]
：圧力損失の二相増倍係数 [-]
：流路計算格子あたりの平均体積 [m3]
：Sherwood 数 [-]
：Nusselt 数 [-]
：Pt 表面代表値
：モデルパラメータ（実測データから決定）
 a
F
R
T

：反応移動係数 [-]
p
：圧力 [Pa]

：反応ガス分圧依存係数 [-]
G
：活性化エネルギー [J/mol]
ilim
：限界電流密度（ガス分圧関数）[A/m2]

：電解質中のプロトン伝導度 [S/m]
w
M
：電解質中の含水量 [mol/mol SO3-]
Multicomponent
：移動モルフラックス [mol/m2/s]
Electrochem. Soc.， 150 (2003), A1589-A1598.
nd
：電解質中の電気浸透係数 [-]
：ファラデー定数 [C/mol]
：ガス定数 [J/mol/K]
：温度 [K]
：触媒利用率（水蒸気活量関数） [-]
： (a) アノード， (c) カソード
：電解質膜
：気相
：液相
：ガス成分
引用文献
1) Wang, C. Y.: Fundamental Models for Fuel Cell
Engineering, Chem. Rev., 104 (2004), 4727-4766.
2) Berning, T. and Djilaji, N.: A 3D, Multiphase,
Model
of
the
Cathode,
J.
3) Meng, H. and Wang, C. Y.: Model of Two-Phase
2
Dw
t
：電解質中の水分拡散係数 [m /s]
kw
EW

u
s

：電解質膜／GDL 界面の水分輸送係数 [m/s]
：粘性係数 [Pa・s]
Mathematical Modeling of PEM Fuel Cells - II.


：多孔質材料の空孔度 [-]
Model Predictions with Liquid Water Transport, J.
Flow and Flooding Dynamics, J. Electrochem. Soc.,
：電解質膜厚 [m]
152 (2005), A1733-A1741.
4) Mazumder, S. and Cole, J. V.: Rigorous 3-D
：電解質の Equivalent Weight [kg/mol SO3-]
Mathematical Modeling of PEM Fuel Cells - I.
3
：密度 [kg/m ]
Model Predictions without Liquid Water Transport, J.
：流速 [m/s]
Electrochem. Soc., 150 (2003), A1503-A1509.
：液相体積率 [m3/m3]
5) Mazumder, S. and Cole, J. V.: Rigorous 3-D
：多孔質材料の屈曲度 [-]
Electrochem. Soc., 150 (2003), A1503-A1509.
2
K
：多孔質材料の透過係数 [m ]
pc
：毛管圧力（液相体積率関数） [Pa]
Ejiri, E.: Macroscopic Modeling and Simulation for
c
w
：ガス濃度 [mol/m3]
Polymer Electrolyte Membrane Fuel Cell, ECS
：分子量 [kg/mol]
Transaction, 3 (2006), 1105-1104.
m p
：電気化学反応による物質生成量 [kg/m3/s]
m v
：相変化による物質生成量 [kg/m3/s]
kv
：相変化速度 [1/s]
6) Yoneda, M., Tago, Y, Suzuki, K., Takimoto, M. and
7) Yoneda,
M.
and
Takimoto,
M.:
Macroscopic
Modeling of Two-Phase Flow Transport Inside
Polymer Electrolyte Fuel Cells, ASME 2010 8th
International Fuel Cell Science, Engineering and
18
Technology Conference, Volume 1 (2010), 497-504.
8) Zhang, F. Y., Yang, X. G. and Wang, C. Y.: Liquid
Water Removal from a Polymer Electrolyte Fuel Cell，
J. Electrochem. Soc., 153 (2006), A225-A232.
9) Yoneda, M., Takimoto, M. and Koshizuka, S.: Effects
of Microstructure of Gas Diffusion Layer on
Two-Phase
Flow
Transport
Properties,
ECS
Transaction, 11 (2007), 1105-1104.
10) Yoneda, M., Takimoto, M., Ejiri, E. and Koshizuka,
S.: Investigation of the Behavior of Liquid Water for
Its Transport Process into the GDL, ECS Transaction,
16 (2008), 787-794.
11) Hirono, M., Suzuma, I., KUrata, T., Tadokoro, R.,
Ejiri,
E.,
Takimoto,
M.
and
Yoneda,
M.:
Measurements of Water Droplet Behavior in PEFC
Channel Using PTV and High Accuracy Pressure
Sensors, ECS Transaction, 25 (2010), 139-145.
12) J. Inukai, K. Miyatake, K. Takada, M. Watanabe, T.
Hyakutake, H. Nishide, Y. Nagumo, M. Watanabe, M.
Aoki, H. Takano, Angew. Chem. Int. Ed. 47 (2008),
2792-2795.
19

Download Report