特許 7304537 相変化を伴う現象のシミュレーション方法及びシミュレーション装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2023-06-29

(45)【発行日】2023-07-07

(54)【発明の名称】相変化を伴う現象のシミュレーション方法及びシミュレーション装置

(51)【国際特許分類】

   H01M 8/04492 20160101AFI20230630BHJP

   H01M 8/04537 20160101ALI20230630BHJP

   H01M 8/04 20160101ALI20230630BHJP

【ＦＩ】

H01M8/04492

H01M8/04537

H01M8/04 Z

【請求項の数】 10

(21)【出願番号】P 2019076485

(22)【出願日】2019-04-12

(65)【公開番号】P2019216085

(43)【公開日】2019-12-19

【審査請求日】2021-11-04

(31)【優先権主張番号】P 2018108649

(32)【優先日】2018-06-06

(33)【優先権主張国・地域又は機関】JP

(73)【特許権者】

【識別番号】314012076

【氏名又は名称】パナソニックＩＰマネジメント株式会社

(74)【代理人】

【識別番号】110000556

【氏名又は名称】弁理士法人有古特許事務所

(72)【発明者】

【氏名】保岡 悠

(72)【発明者】

【氏名】山本 恵一

(72)【発明者】

【氏名】藤田 悠介

(72)【発明者】

【氏名】松下 史弥

【審査官】藤森 一真

(56)【参考文献】

【文献】特開２００９－２４５８７４（ＪＰ，Ａ）

【文献】特開２００９－１９３６７２（ＪＰ，Ａ）

【文献】特開２０１７－１３０３０６（ＪＰ，Ａ）

【文献】特開２００９－２１１９４０（ＪＰ，Ａ）

【文献】特開２００５－１３５８１４（ＪＰ，Ａ）

【文献】特開２０１０－１６０９９３（ＪＰ，Ａ）

【文献】L. Matamoros, D. Bruggemann，Simulation of water and heat management in proton exchange membrane fuel cells，Journal of Power Sources，Elsevier，2006年05月18日，Vol.161，p.203-213

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｈ０１Ｍ ８／０４ － ８／０６６８

(57)【特許請求の範囲】

【請求項1】

多孔質体における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量に基づいて、前記多孔質体を含む燃料電池の液水分布を算出し、
そこでは、前記多孔質体において液水以外の空隙に拡散するガスの実効的な相互拡散係数を、ボルツマン方程式の考え方に基づくシミュレーション手法を用いて算出し、
前記ガスの実効的な相互拡散係数に基づいて、前記燃料電池における前記ガスの分布を算出し、前記燃料電池における前記液水分布及び前記ガスの分布に基づいて、前記燃料電池における電位の分布を算出する、
相変化を伴う現象のシミュレーション方法であって、
前記液水分布の算出に式４を用い、

【数1】

前記ガスの実効的な相互拡散係数の算出に式５で示されるＳｔｅｆａｎ－Ｍａｘｗｅｌｌ方程式を用いる、

【数2】

相変化を伴う現象のシミュレーション方法。
（式４中、εは細孔に対する空隙の割合（空隙率）であり、ρ _ｌ は液水の密度であり、ｓは細孔内の液水の充填率であり、ｔは時間であり、Ｋは液水の透過率であり、Ｐｃは毛管圧であり、μ _ｌ は液水の粘性係数であり、ｒ _ｗ は液水の単位体積、単位時間当たりの相変化量である。）
（式５中、Ｘ _ｉ は第１ガス粒子ｉのモル分率であり、Ｘ _ｊ は第２ガス粒子ｊのモル分率であり、Ｖ _ｉ は第１ガス粒子ｉの速度であり、Ｖ _ｊ は第２ガス粒子ｊの速度であり、Ｄ _ｉｊ ^ｅｆｆ は、第１ガス粒子ｉと第２ガス粒子ｊとの実効的な相互拡散係数（ガスの実効的な相互拡散係数）である。）

【請求項2】

前記相変化量の算出前に、
前記多孔質体の形状及び材質に関する情報を含む構造情報を取得し、
前記多孔質体の構造情報に基づき、前記多孔質体における単位体積当たりの気液界面の面積を算出する、
請求項１記載の相変化を伴う現象のシミュレーション方法。

【請求項3】

前記相変化量の算出前に、前記多孔質体における単位体積当たりの気液界面の面積を取得する、
請求項１に記載の相変化を伴う現象のシミュレーション方法。

【請求項4】

前記相変化量の算出前に、分子動力学法を用いて前記気液界面の相変化速度を算出する、
請求項１～３のいずれか一項に記載の相変化を伴う現象のシミュレーション方法。

【請求項5】

前記相変化量の算出前に、前記気液界面の相変化速度を取得する、
請求項１～３のいずれか一項に記載の相変化を伴う現象のシミュレーション方法。

【請求項6】

演算処理部を備え、
前記演算処理部は、多孔質体における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量に基づいて、前記多孔質体を含む燃料電池の液水分布を算出し、
そこでは、前記多孔質体において液水以外の空隙に拡散するガスの実効的な相互拡散係数を、ボルツマン方程式の考え方に基づくシミュレーション手法を用いて算出し、
前記ガスの実効的な相互拡散係数に基づいて、前記燃料電池における前記ガスの分布を算出し、前記燃料電池における前記液水分布及び前記ガスの分布に基づいて、前記燃料電池における電位の分布を算出するように構成されている、
相変化を伴う現象のシミュレーション装置であって、
前記演算処理部は、前記液水分布の算出に式４を用い、

【数3】

前記ガスの実効的な相互拡散係数の算出に式５で示されるＳｔｅｆａｎ－Ｍａｘｗｅｌｌ方程式を用いるように構成されている、

【数4】

相変化を伴う現象のシミュレーション装置。
（式４中、εは細孔に対する空隙の割合（空隙率）であり、ρ _ｌ は液水の密度であり、ｓは細孔内の液水の充填率であり、ｔは時間であり、Ｋは液水の透過率であり、Ｐｃは毛管圧であり、μ _ｌ は液水の粘性係数であり、ｒ _ｗ は液水の単位体積、単位時間当たりの相変化量である。）
（式５中、Ｘ _ｉ は第１ガス粒子ｉのモル分率であり、Ｘ _ｊ は第２ガス粒子ｊのモル分率であり、Ｖ _ｉ は第１ガス粒子ｉの速度であり、Ｖ _ｊ は第２ガス粒子ｊの速度であり、Ｄ _ｉｊ ^ｅｆｆ は、第１ガス粒子ｉと第２ガス粒子ｊとの実効的な相互拡散係数（ガスの実効的な相互拡散係数）である。）

【請求項7】

前記演算処理部は、
前記相変化量の算出前に、
前記多孔質体の形状及び材質に関する情報を含む構造情報を取得し、
前記多孔質体の構造情報に基づき、前記多孔質体における単位体積当たりの気液界面の面積を算出するように構成されている、
請求項６記載の相変化を伴う現象のシミュレーション装置。

【請求項8】

前記演算処理部は、前記相変化量の算出前に、前記多孔質体における単位体積当たりの気液界面の面積を取得するように構成されている、
請求項６に記載の相変化を伴う現象のシミュレーション装置。

【請求項9】

前記演算処理部は、前記相変化量の算出前に、分子動力学法を用いて前記気液界面の相変化速度を算出するように構成されている、
請求項６～８のいずれか一項に記載の相変化を伴う現象のシミュレーション装置。

【請求項10】

前記演算処理部は、前記相変化量の算出前に、前記気液界面の相変化速度を取得するように構成されている、
請求項６～８のいずれか一項に記載の相変化を伴う現象のシミュレーション装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、相変化を伴う現象のシミュレーション方法に関する。

【背景技術】

【0002】

多孔質体を含む燃料電池では、多孔質体における液体の水（液水）の輸送により、その性能が影響を受ける。このため、例えば、特許文献１の燃料電池のシミュレーション方法では、触媒層の構造を算出し、水の輸送等に基づき触媒層の性能を算出し、触媒層の性能を基に燃料電池の発電性能を算出することにより、電解質膜と電極とを有する膜－電極を備える燃料電池の発電性能をシミュレーションしている。

【0003】

また、多孔質体における液水の輸送は、液水に加わる圧力勾配に加え、蒸発及び凝縮の相変化にも影響を受ける。この液水の蒸発量及び凝縮量の相変化量について、例えば、非特許文献１の方法によれば、多孔質体における液水の充填率（飽和度）に比例するように決定されている。

【先行技術文献】

【特許文献】

【0004】

【文献】特開２０１５－１３５７９０号公報

【非特許文献】

【0005】

【文献】Ｅｎｅｒｇｙ Ｖｏｌ． ７３ ｐｐ６１８－６３４ （２０１４）

【発明の概要】

【発明が解決しようとする課題】

【0006】

本発明は、高精度な相変化を伴う現象のシミュレーション方法の構築を課題としている。

【課題を解決するための手段】

【0007】

本発明のある態様に係る相変化を伴う現象のシミュレーション方法では、多孔質体における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量を算出する。

【発明の効果】

【0008】

本発明は、高精度な相変化を伴う現象のシミュレーション方法を提供することができるという効果を奏する。

【0009】

本発明の上記目的、他の目的、特徴、及び利点は、添付図面参照の下、以下の好適な実施態様の詳細な説明から明らかにされる。

【図面の簡単な説明】

【0010】

【図1】本発明の実施の形態１に係る相変化を伴う現象のシミュレーション装置を示す機能ブロック図である。

【図2】相変化を伴う現象のシミュレーション方法に用いられる多孔質体を模式的に示した断面図である。

【図3】気液界面の相変化速度に関する解析領域を模式的に示す斜視図である。

【図4】気体分子の数の経時変化を示すグラフである。

【図5】その他の変形例に係るシミュレーション方法の一例を示すフローチャートである。

【図6】燃料電池を模式的に示す斜視図である。

【図7】本発明の実施の形態２に係る相変化を伴う現象のシミュレーション方法の一例を示すフローチャートである。

【図8】ガスの相互拡散係数に関する解析領域を模式的に示す斜視図である。

【図9】本発明の実施の形態３に係る相変化を伴う現象のシミュレーション方法の一例を示すフローチャートである。

【図10】実施の形態１に係るシミュレーション方法の一例を示すフローチャートである。

【図11】変形例１に係るシミュレーション方法の一例を示すフローチャートである。

【図12】変形例２に係るシミュレーション方法の一例を示すフローチャートである。

【図13】変形例３に係るシミュレーション方法の一例を示すフローチャートである。

【図14】変形例４に係るシミュレーション方法の一例を示すフローチャートである。

【発明を実施するための形態】

【0011】

（本発明の基礎となった知見）
本発明者らは、多孔質体における液水に関するシミュレーションの精度向上について検討した。液水の相変化は、液水の気液界面を介して水分子が移動することにより生じる。しかしながら、非特許文献１の方法では、液水の相変化量は、多孔質体における液水の充填率に基づいており、多孔質体における気液界面の面積が液水の相変化量に反映されていない。そこで、本発明者等は、多孔質体における気液界面の面積を考慮して液水の相変化量を求めることにより、シミュレーションの精度向上が図られることを見出した。

【0012】

第１の発明に係る相変化を伴う現象のシミュレーション方法では、多孔質体における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量を算出する。

【0013】

この方法によれば、液水の相変化量（蒸発量及び凝縮量）は、多孔質体における単位体積当たりの気液界面の面積に基づいて算出される。このため、液水の相変化量について、相変化が生じる気体と液水の気液界面を考慮した高精度なシミュレーションを行うことができる。

【0014】

第２の発明に係る相変化を伴う現象のシミュレーション方法では、第１の発明において、前記相変化量の算出前に、前記多孔質体の構造情報を取得し、前記多孔質体の構造情報に基づき、前記多孔質体における単位体積当たりの気液界面の面積を算出する。この方法によれば、多孔質体の構造を考慮した液水の相変化量のシミュレーションを行うことができる。

【0015】

第３の発明に係る相変化を伴う現象のシミュレーション方法では、第１の発明において、前記相変化量の算出前に、前記多孔質体における単位体積当たりの気液界面の面積を取得する。この方法によれば、取得した気液界面の面積を考慮した液水の相変化量を高精度にシミュレーションすることができる。

【0016】

第４の発明に係る相変化を伴う現象のシミュレーション方法では、第１～３のいずれかの発明において、前記相変化量の算出前に、分子動力学法を用いて前記気液界面の相変化速度を算出する。この方法によれば、単位体積当たりの単位面積当たりに蒸発する液水量を決定する相変化速度係数を決定することができる。この方法によれば、分子動力学法を用いて算出した相変化速度を考慮した液水の相変化量を高精度にシミュレーションすることができる。また、このシミュレーション方法においては、分子動力学法を用いることにより、実験等に比べて容易に気液界面の相変化速度を取得することができる。

【0017】

第５の発明に係る相変化を伴う現象のシミュレーション方法では、第１～３のいずれかの発明において、前記相変化量の算出前に、前記気液界面の相変化速度を取得する。この方法によれば、取得した相変化速度を考慮した液水の相変化量を高精度にシミュレーションすることができる。

【0018】

第６の発明に係る相変化を伴う現象のシミュレーション方法では、第１～５のいずれかの発明において、前記多孔質体の単位体積、単位時間当たりの相変化量に基づいて、前記多孔質体を含む燃料電池の液水分布を算出する。この方法によれば、燃料電池における液水分布について、液水の気液界面を考慮した高精度なシミュレーションを行うことができる。

【0019】

第７の発明に係る相変化を伴う現象のシミュレーション方法では、第６の発明において、前記多孔質体において液水以外の空隙に拡散するガスの実効的な相互拡散係数を、ボルツマン方程式の考え方に基づくシミュレーション手法を用いて算出し、前記ガスの実効的な相互拡散係数に基づいて、前記燃料電池における前記ガスの分布を算出し、前記燃料電池における前記液水分布及び前記ガスの分布に基づいて、前記燃料電池における電位の分布を算出する。

【0020】

この方法によれば、燃料電池の電位分布（発電性能予測）のシミュレーションを高精度に行うことができる。なお、ボルツマン方程式の考え方に基づくシミュレーション手法とは、格子ボルツマン法及び格子ガス法のことを意味している。

【0021】

以下、本発明の実施形態を、図面を参照しながら具体的に説明する。なお、以下では全ての図面を通じて同一又は相当する要素には同一の参照符号を付して、その重複する説明を省略する。

【0022】

（実施の形態１）
＜シミュレーション装置の構成＞
本発明の実施の形態１に係る相変化を伴う現象のシミュレーション装置１０は、図１に示すように、多孔質体における液体の水（液水）の相変化量等を予測する装置である。この多孔質体には、燃料電池セルに用いられるガス拡散層及び触媒層等が例示される。このため、シミュレーション装置１０は、例えば、燃料電池におけるシミュレーション装置として用いられるが、これ以外のシミュレーション装置にも用いることができる。

【0023】

液水の相変化は、水の液体から気体への変化（蒸発）、及び、水の気体から液体への変化（凝縮）である。また、液水の相変化量は、蒸発する水の量（蒸発量）、又は、凝縮する水の量（凝縮量）である。

【0024】

シミュレーション装置１０は、演算処理部１１及び記憶部１２を有しており、例えば、コンピュータによって構成されている。例えば、演算処理部１１としてコンピュータのＣＰＵ等が例示され、記憶部１２として演算処理部１１によりアクセス可能なメモリ等が例示される。

【0025】

記憶部１２には、多孔質体における液水の相変化等をシミュレーションするプログラム及びこれに必要な情報が記憶されている。なお、このプログラム及びこれに必要な情報は、コンピュータに内蔵されたメモリに記憶されたものに限定されず、他の記憶媒体に記憶されたものであってもよいし、入力装置によって入力されたものであってもよいし、ネットワークを介して受信されたものであってもよい。

【0026】

演算処理部１１は、記憶部１２に記憶されたプログラム等のソフトウェアを読み出して実行することにより、多孔質体における液水の相変化等をシミュレーションする。シミュレーション装置１０は、それぞれ、単独の装置で構成されていても構わないし、互いに協働する複数の装置で構成されていても構わない。

【0027】

＜多孔質体の構成＞
シミュレーションの対象とする図２の多孔質体２０について説明する。多孔質体２０には、壁部分２１及び多数の細孔２２が設けられている。壁部分２１は、例えば、樹脂及びカーボン担体等の有機物、ガラス等の無機物、並びに、その混合物等により形成されている。

【0028】

細孔２２は、多孔質体２０において周囲が壁部分２１により取り囲まれていた空間である。細孔２２に液体状の水（液水２４）が分布している場合、細孔２２から液水２４を除いた空間（空隙２３）が、ガスの粒子が移動可能な空間となる。このため、空隙２３は、多孔質体２０において壁部分２１及び液水２４以外の領域である。例えば、燃料電池の場合、ガスの粒子としては、水素、酸素及び窒素等が例示される。

【0029】

＜相変化を伴う現象のシミュレーション方法＞
相変化を伴う現象のシミュレーション方法は、例えば、図１のシミュレーション装置１０によって図１０に示すフローチャートに沿って実行される。ここでは、図２に示す多孔質体２０の単位体積及び単位時間当たりの液水２４の相変化量についてシミュレーションする。

【0030】

シミュレーション装置１０の演算処理部１１は、多孔質体２０における単位体積当たりの気液界面の面積、及び、気液界面の相変化速度から、多孔質体２０の単位体積、単位時間当たりの相変化量を算出する（ステップＳ２０）。

【0031】

この液水２４の相変化は、多孔質体２０の細孔２２において空隙２３に接する液水２４の気液界面（液水２４の表面）において生じる。多孔質体２０における単位体積当たりの気液界面の面積（比表面積Ｓ_ｆ）は、単位体積の多孔質体２０において細孔２２の空隙２３に接している液水２４の気液界面の面積である。また、気液界面の相変化速度（相変化速度係数ｋ）は、気液界面において液水２４が蒸発する速度、又は、気液界面において液水２４に凝縮する速度である。

【0032】

続いて、演算処理部１１は、相変化速度係数ｋ及び比表面積Ｓ_ｆに基づき、下記式３から多孔質体２０の単位体積及び単位時間当たりの液水２４の相変化量（多孔質体２０の単位体積当たりの気液界面の相変化速度（単位体積、単位時間当たりの相変化量ｒ_ｗ））を算出する。

【0033】

この式３において、ｋは、相変化速度係数である。液水２４の水蒸気分圧ｐ_ｗｖ、液水２４の飽和蒸気圧ｐ_ｓａｔ、温度Ｔは、求めたいシミュレーション条件により定められる。

【0034】

【数1】

【0035】

このような相変化を伴う現象のシミュレーション方法では、液水２４の比表面積Ｓ_ｆ及び相変化速度係数ｋから、単位体積、単位時間当たりの相変化量ｒ_ｗを算出している。このように、液水２４の比表面積Ｓ_ｆによって、多孔質体２０において相変化が生じる液水２４の気液界面の面積（表面積）を考慮して、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗを高精度にシミュレーションすることができる。

【0036】

＜変形例１＞
変形例１に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、多孔質体２０の構造情報を取得し、多孔質体２０の構造情報に基づき、多孔質体２０における単位体積当たりの気液界面の面積を算出する。このシミュレーション方法は、図１のシミュレーション装置１０によって図１１に示すフローチャートに沿って実行される。

【0037】

具体的には、演算処理部１１は、図１１のシミュレーション方法において対象とする多孔質体２０の構造情報を取得する（ステップＳ２１）。この多孔質体２０の構造情報は、多孔質体２０の形状及び材質に関する情報を含む。

【0038】

例えば、多孔質体２０の壁部分２１の材質情報としては、壁部分２１に対する液水２４の接触角θが例示され、これは実験等より予め定めされうる。また、多孔質体２０内の液水２４の情報として、表面張力γがあり、これは実験等より予め定めされている。

【0039】

壁部分２１の形状情報は、多孔質体２０の画像情報等により得られる。例えば、ＦＩＢ－ＳＥＭ（Ｆｏｃｕｓｅｄ Ｉｏｎ Ｂｅａｍ－Ｓｃａｎｎｉｎｇ Ｅｌｅｃｔｒｏｎ Ｍｉｃｒｏｓｃｏｐｅ）を用いて、多孔質体２０について所定の方向において異なる位置で複数の断面画像を連続的に撮影する。

【0040】

この複数の断面画像を積層した画像（積層画像）を画像処理によって二値化の作業を行う。これにより、細孔面２０ａの位置情報を含むシミュレーション用の計算メッシュが、多孔質体２０の形状情報として作成される。ここで、積層画像に基づいた多孔質体２０の形状情報の取得には、Ｍａｔｈ２Ｍａｒｋｅｔ社製のＧｅｏＤｉｃｔ（ドイツの登録商標）等の解析ソフトウェアを用いてもよい。

【0041】

そして、演算処理部１１は、ステップＳ２１において情報が取得された多孔質体２０における液水２４の分布（液水分布）を取得する（ステップＳ２２）。液水分布は、多孔質体２０の細孔２２に配置される液水２４の位置情報（座標等）であって、以下のとおり、式２及びポアモルフォロジー法による液水２４を詰め込むシミュレーションを行うことにより求められる。

【0042】

【数2】

【0043】

例えば、これは、Ｍａｔｈ２Ｍａｒｋｅｔ社製の解析ソフトＧｅｏＤｉｃｔ等のソフトウェアに用いることにより算出することが可能である。具体的には、まず、多孔質体２０の壁部分２１の材質情報、及び、圧力（毛管圧）により、多孔質体２０内に液水２４が存在可能な空隙の半径をヤング・ラプラスの式（式２）により決定する。式２において、ｐ_ｃは毛管圧［Ｐａ］である。γは液水２４の表面張力［Ｎ／ｍ］であり、θは壁部分２１に対する液水２４の接触角［°］であって、これらは壁部分２１の材質情報から求められる。ｒは、空隙の半径［ｍ］である。

【0044】

多孔質体２０における壁部分２１に対する液水２４の接触角θが９０°を超える（θ＞９０°）場合、壁部分２１は撥水性である。このため、液水２４は、毛管圧Ｐ_ｃによって細孔２２から押し出される力を受ける。毛管圧Ｐ_ｃから式２に基づいて液水２４が浸入できる空隙２３の最小の半径ｒを求める。そして、多孔質体２０の所定面から液水２４を順次、詰め込むシミュレーションを、空隙２３の最小の半径ｒを基にポアモルフォロジー法により実行する。これにより、多孔質体２０に配置される液水２４の位置情報が液水分布として求められる。

【0045】

一方、液水２４の接触角θが９０°を下回る（θ＜９０°）場合、多孔質体２０の壁部分２１は親水性である。このため、液水２４は、毛管圧Ｐ_ｃによって細孔２２に取り込まれる力を受ける。毛管圧Ｐ_ｃから式２に基づいて気体により液水２４が排出される空隙２３の最大の半径ｒを求める。そして、上記多孔質体２０の形状情報に基づく多孔質体２０の所定面から液水２４が記載により順次、押し出されるシミュレーションを、空隙２３の最大の半径ｒを基に上記ポアモルフォロジー法により実行する。これにより、多孔質体２０に配置される液水２４の位置情報が液水分布として求められる。

【0046】

ここで、毛管圧Ｐ_ｃを変化させることにより、多孔質体２０における液水２４の充填率が変化する。このため、式２において毛管圧Ｐ_ｃを変えながら、シミュレーションを実行することにより、液水分布が充填率の関数として得られる。なお、液水２４の充填率は０から１の範囲に設定される。液水２４の充填率が０は、細孔２２に液水２４が配置されていない状態を示し、液水２４の充填率が１は、細孔２２に液水２４が満たされている状態を示す。

【0047】

例えば、Ｍａｔｈ２Ｍａｒｋｅｔ社製の解析ソフトＧｅｏＤｉｃｔを用いることにより、多孔質体２０の構造情報及び液水２４の充填率から液水分布が決定される。

【0048】

続いて、演算処理部１１は、ステップＳ２２の液水分布に基づき、多孔質体２０の単位体積当たりの気液界面の面積（比表面積Ｓ_ｆ）を算出する（ステップＳ２３）。

【0049】

ここで、Ｓ２１における壁部分２１の形状情報から、細孔２２の位置が求められる。また、Ｓ２２における液水分布から、細孔２２における液水２４の位置が求められる。このため、これらの位置から、多孔質体２０における空隙２３に接する液水２４の気液界面の面積（表面積）を算出する。

【0050】

そして、気液界面の面積を多孔質体２０の体積により割ることにより、液水２４の比表面積Ｓ_ｆが算出される。ここで、液水分布が充填率の関数で表されている場合、これに基づく液水２４の比表面積Ｓ_ｆも充填率の関数で定められる。例えば、液水２４の比表面積Ｓ_ｆは、Ｍａｔｈ２Ｍａｒｋｅｔ社製のＧｅｏＤｉｃｔ等のソフトウェアを用いることにより得られる。

【0051】

続いて、演算処理部１１は、気液界面の相変化速度係数ｋ、及び、ステップＳ２３の液水２４の比表面積Ｓ_ｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒ_ｗを算出する（ステップＳ２０）。

【0052】

このように、多孔質体２０の構造情報に基づいて算出した液水２４の比表面積Ｓ_ｆを用いて、多孔質体２０の単位体積、単位時間当たりの相変化量ｒ_ｗを算出している。これにより、多孔質体２０の構造を考慮した、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗを高精度にシミュレーションすることができる。

【0053】

＜変形例２＞
変形例２に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、多孔質体２０における単位体積当たりの気液界面の面積を取得する。このシミュレーション方法は、図１１に示すフローチャートにおけるステップＳ２１～Ｓ２３の処理に代えて、図１２に示すフローチャートにおけるステップＳ２４の処理が実行される。

【0054】

具体的には、演算処理部１１は、図１２のシミュレーション方法において対象とする多孔質体２０の単位体積当たりの気液界面の面積（比表面積Ｓ_ｆ）を取得する（ステップＳ２４）。この比表面積Ｓ_ｆは、シミュレーション装置１０に対する入力値であってもよく、例えば、記憶部１２に記憶されたものに限定されず、他の記憶媒体に記憶されたものであってもよいし、入力装置によってシミュレーション装置１０に入力されたものであってもよいし、ネットワークを介して受信されたものであってもよい。

【0055】

そして、演算処理部１１は、気液界面の相変化速度係数ｋ、及び、ステップＳ２４の液水２４の比表面積Ｓ_ｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒ_ｗを算出する（ステップＳ２０）。

【0056】

このように、取得した液水２４の比表面積Ｓ_ｆを用いて、多孔質体２０の単位体積、単位時間当たりの相変化量ｒ_ｗを算出する。これにより、取得した比表面積Ｓ_ｆに応じた、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗを高精度にシミュレーションすることができる。

【0057】

＜変形例３＞
変形例３に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、分子動力学法を用いて気液界面の相変化速度を算出する。このシミュレーション方法は、図１のシミュレーション装置１０によって図１３に示すフローチャートに沿って実行される。

【0058】

具体的には、演算処理部１１は、気液界面の相変化速度（相変化速度係数ｋ）を算出する（ステップＳ２５）。この相変化の速度（相変化速度）は、例えば、図３に示す解析領域３０において分子動力学法によるシミュレーションを行うことによって求められる。

【0059】

解析領域３０は、例えば、直方体形状の空間である。この解析領域３０に、膜のように薄い直方体形状の液水２４（液膜３１）が、配置されている。液膜３１は、解析領域３０の長手方向に対して直交する方向に拡がり、解析領域３０を長手方向に２つの部分（第１空間３２、第２空間３３）に区分けするように配置されている。このため、液膜３１は、第１空間３２に接する一方面（第１主面３１ａ）、及び、第２空間３３に接する他方面（第２主面３１ｂ）を有している。第１主面３１ａ及び第２主面３１ｂから成る液水２４の気液界面において液膜３１の水分子が相変化する。

【0060】

解析領域３０において平衡状態の液膜３１に熱エネルギーを加える。これにより、液膜３１の水分子が液水２４の気液界面から蒸発し、気体分子となる。この気体分子の数をカウントする。この気体分子の数の経時変化が図４のグラフに求められる。

【0061】

図４のグラフに示すように、時刻０から時刻ａまでは気体分子の数の経時変化が所定値以下であり、液膜３１が平衡状態である。そして、時刻ａで所定量のエネルギーを液膜３１に加える。この時刻ａにおいて気体分子の数の経時変化が増加し、水分子の蒸発が開始する。そして、時刻ｂにおいて気体分子の数の経時変化が減少し、水分子の蒸発が終了して、液膜３１が平衡状態となる。この時刻ａから時刻ｂまでの気体分子の数の経時変化から、気液界面の相変化速度ｒ_０が求められる。

【0062】

この気液界面の相変化速度ｒ_０及びシミュレーション条件を下記式１に代入することにより、気液界面の単位面積当たりの液水２４の相変化速度（相変化速度係数ｋ）が算出される。この式１において、Ｓ_０は、単位体積当たりの液水２４の気液界面の面積（液水２４の表面積）であって、液膜３１の第１主面３１ａ及び第２主面３１ｂの合計面積を解析領域３０の体積で割った値である。ｐｗｖは、液水２４の水蒸気分圧であり、ｐ_ｓａｔは、液水２４の飽和蒸気圧であり、Ｔは、温度であって、これらはシミュレーション条件により定められる。Ｒは、気体定数であり、ＭＨ_２Ｏは、液水２４の分子量である。

【0063】

【数3】

【0064】

そして、演算処理部１１は、ステップＳ２５の気液界面の相変化速度係数ｋ、及び、液水２４の比表面積Ｓ_ｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒ_ｗを算出する（ステップＳ２０）。

【0065】

このように、算出した気液界面の相変化速度を用いて、多孔質体２０の単位体積、単位時間当たりの相変化量ｒ_ｗを算出している。これにより、算出条件に応じた、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗを高精度にシミュレーションすることができる。

【0066】

また、このシミュレーション方法においては、分子動力学法によるシミュレーションによって気液界面の相変化速度を取得可能である。この相変化速度から気液界面の相変化速度係数ｋを実験等に比べて容易に取得することができる。

【0067】

なお、ステップＳ２０において、液水２４の比表面積Ｓ_ｆは、例えば、変形例１のようにシミュレーション装置１０によって算出されてもよいし、変形例２のようにシミュレーション装置１０によって取得されてもよい。

【0068】

例えば、気液界面の相変化速度係数ｋ、及び、液水２４の比表面積Ｓ_ｆがシミュレーション装置１０によって算出される場合、相変化を伴う現象のシミュレーション方法は、図５に示すフローチャートに沿って実行される。この図５のステップＳ１は図１３のステップＳ２５に相当し、図５のステップＳ２～Ｓ４は図１１のステップＳ２１～Ｓ２３にそれぞれ相当し、図５のステップＳ５は図１１のステップＳ２０に相当する。

【0069】

＜変形例４＞
変形例４に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、気液界面の相変化速度を取得する。このシミュレーション方法は、図１３に示すフローチャートにおけるステップＳ２５の処理に代えて、図１４に示すフローチャートにおけるステップＳ２６の処理が実行される。

【0070】

具体的には、演算処理部１１は、気液界面の相変化速度（相変化速度係数ｋ）を取得する（ステップＳ２５）。この相変化速度係数ｋは、シミュレーション装置１０に対する入力値であってもよく、例えば、記憶部１２に記憶されたものに限定されず、他の記憶媒体に記憶されたものであってもよいし、入力装置によって入力されたものであってもよいし、ネットワークを介して受信されたものであってもよい。

【0071】

そして、演算処理部１１は、ステップＳ２６の気液界面の相変化速度係数ｋ、及び、液水２４の比表面積Ｓ_ｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒ_ｗを算出する（ステップＳ２０）。

【0072】

このように、取得した相変化速度係数ｋに応じた、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗを高精度にシミュレーションすることができる。

【0073】

なお、ステップＳ２０において、液水２４の比表面積Ｓ_ｆは、例えば、変形例１のようにシミュレーション装置１０によって算出されてもよいし、変形例２のようにシミュレーション装置１０によって取得されてもよい。

【0074】

（実施の形態２）
本発明の実施の形態２に係るシミュレーション装置１０は、図６に示す燃料電池５０における液水２４の分布（液水分布）をその発電性能として予測する装置である。

【0075】

燃料電池５０は、例えば、電解質膜５１、これを互いの間に挟む一対の触媒層５２、これを互いの間に挟む一対のガス拡散層５３、及び、これを互いの間に挟むガス流路５４により構成されている。このうち、触媒層５２及びガス拡散層５３の少なくともいずれか１つの部材が多孔質体２０により形成されている。

【0076】

シミュレーション装置１０による相変化を伴う現象のシミュレーション方法では、燃料電池５０内における液水２４の充填率を算出する。

【0077】

＜燃料電池における液水分布の決定＞
多孔質体２０の単位体積当たりの相変化速度（単位体積、単位時間当たりの相変化量ｒ_ｗ）に基づいて、燃料電池５０における液水分布を下記式４により算出することができる。ここで、燃料電池５０における条件を式４に代入する。

【0078】

式４において、εは細孔２２に対する空隙２３の割合（空隙率）であり、ρ_ｌは液水２４の密度であり、ｓは細孔２２内の液水２４の充填率であり、ｔは時間［ｓｅｃ］である。また、Ｋは液水２４の透過率であり、ｒは液水２４の充填率の次数である。Ｐｃは毛管圧であり、上記式１により求められる。μ_ｌは液水２４の粘性係数である。ｒ_ｗは液水２４の単位体積、単位時間当たりの相変化量であり、上記式３から求められる。

【0079】

【数4】

【0080】

このように、燃料電池５０内での液水２４の充填率ｓの分布が式４により求められる。これは、例えば、ＦＬＵＥＮＴ（登録商標）等のソフトウェアに導入することにより算出される。

【0081】

＜相変化を伴う現象のシミュレーション方法＞
相変化を伴う現象のシミュレーション方法は、例えば、シミュレーション装置１０によって図７に示すフローチャートに沿って実行される。図７に示すフローチャートでは、図１０に示すステップＳ２０の処理に加えて、ステップＳ６の処理をさらに実行する。なお、図７のステップＳ２０の処理の前に、図１１のステップＳ２１～Ｓ２３の処理、又は、図１２のステップＳ２４の処理を行ってもよい。また、図７のステップＳ２０の処理の前に、図１３のステップＳ２５の処理、又は、図１４のステップＳ２６の処理を行ってもよい。

【0082】

シミュレーション装置１０の演算処理部１１は、多孔質体２０の単位体積当たりの相変化速度（単位体積、単位時間当たりの相変化量ｒ_ｗ）を算出する（ステップＳ２０）。ここでは、ステップＳ２０の処理により液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗが求められる。

【0083】

また、演算処理部１１は、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗに基づいて、燃料電池５０における液水２４の分布（液水分布）を算出する（ステップＳ６）。この燃料電池５０における液水分布は、上述の通り、式４により求められる。

【0084】

この相変化を伴う現象のシミュレーション方法では、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗに基づいて燃料電池５０における液水分布を算出している。これにより、液水２４の比表面積Ｓ_ｆに基づいた液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗによって、多孔質体２０において相変化が生じる液水２４の気液界面の面積を考慮して、燃料電池５０における液水分布を高精度にシミュレーションすることができる。

【0085】

また、式４における液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗに、多孔質体２０における液水２４の輸送が依存する。このため、多孔質体２０から成る部材を含む燃料電池５０において、液水２４の輸送を考慮した高精度のシミュレーションを行うことができる。

【0086】

（実施の形態３）
本発明の実施の形態３に係るシミュレーション装置１０は、図６に示す燃料電池５０の電位の分布をその発電性能として予測する装置である。燃料電池５０の電位分布は、燃料電池５０における液水２４の分布（液水分布）及びガスの分布に基づいて求められる。燃料電池５０におけるガスの分布は、ガスの実効的な相互拡散係数に基づいて算出される。

【0087】

＜ガスの実効的な相互拡散係数の決定＞
ガスの粒子が細孔２２における空隙２３において移動しながら、ガスが拡散する。このガスの拡散について、ボルツマン方程式の考え方に基づくシミュレーション手法を用いたガス拡散シミュレーションでは、図８に示す流路構造の解析領域４０が用いられる。

【0088】

解析領域４０は、多孔質体２０、第１流路４１及び第２流路４２を有している。第１流路４１及び第２流路４２は、例えば、断面が矩形状であって、第１方向に長く延びており、互いに間隔を空けて互いに平行に設けられている。多孔質体２０は、例えば、直方体形状であって、第１流路４１と第２流路４２との間に配置されており、第１方向に直交する第２方向に延びている。

【0089】

第１方向における多孔質体２０の寸法は第１流路４１及び第２流路４２の各寸法よりも小さい。第２方向における多孔質体２０の寸法は第１流路４１及び第２流路４２の各寸法よりも大きい。第１方向及び第２方向に直交する第３方向における孔質体の寸法は第１流路４１及び第２流路４２の各寸法と等しい。

【0090】

このような解析領域４０において、多孔質体２０は第１流路４１及び第２流路４２に接続されており、多孔質体２０の細孔２２は第１流路４１及び第２流路４２とそれぞれ連通している。細孔２２は第１流路４１側の面と第２流路４２側の面との間に長く延びており、第１流路４１と第２流路４２とは細孔２２を介して互いに連通している。

【0091】

ここで、液水２４が多孔質体２０に充填され、細孔２２に形成された空隙２３と第１流路４１及び第２流路４２とが連通する。この第１流路４１及び第２流路４２には互いに異なる種類のガス（第１ガス、第２ガス）を流入する。この流入するガスの流速は、実際の燃料電池の触媒層の中を流れるガスの平均的な流速に設定する。

【0092】

第１ガスの粒子（第１ガス粒子ｉ）は第１流路４１の入口から流入し、その一部が多孔質体２０の空隙２３に流れ込んで拡散し、残りは第１流路４１の出口から流出する。第２ガスの粒子（第２ガス粒子ｊ）は第２流路４２の入口から流入し、その一部が多孔質体２０の空隙２３に流れ込んで拡散し、残りは第２流路４２の出口から流出する。

【0093】

これにより、第１流路４１の出口における第１ガスのモル分率、及び、第２流路４２の出口における第２ガスのモル分率は、空隙２３における第１ガス及び第２ガスの拡散状態（ガスの相互拡散）に依存する。このため、上記格子ボルツマン法のガス拡散シミュレーションによる各流路の出口における第１ガスのモル分率及び第２ガスのモル分率が、下記式５で示されるＳｔｅｆａｎ－Ｍａｘｗｅｌｌ方程式のモル分率に一致するように、実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆを調整することによって、実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆが求められる。例えば、これは、ＦＬＵＥＮＴ等の熱流体解析ソフトウェアを利用することにより算出することが可能である。

【0094】

【数5】

【0095】

式５において、Ｘ_ｉは第１ガス粒子ｉのモル分率であり、Ｘ _ｊ は第２ガス粒子ｊのモル分率である。Ｖ_ｉは第１ガス粒子ｉの速度であり、Ｖ_ｊは第２ガス粒子ｊの速度である。Ｄ_ｉｊ ^ｅｆｆは、第１ガス粒子ｉと第２ガス粒子ｊとの実効的な相互拡散係数（ガスの実効的な相互拡散係数）である。ガスの実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆは、多孔質体２０の形状及び多孔質体２０における液水分布を考慮したガスの相互拡散係数である。

【0096】

この多孔質体２０における液水分布は、多孔質体２０における液水２４の充填率に依存する。よって、ガスの実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆは、充填率を変数とする関数として得られる。

【0097】

＜燃料電池の電位分布の決定＞
燃料電池５０の電位は、下記式６、式７及び式８で表される電子及びプロトンに関する方程式により算出される。また、例えば、ＦＬＵＥＮＴ等のソフトウェアに用いることにより算出されうる。

【0098】

【数6】

【0099】

【数7】

【0100】

式６において、σ_ｅ ^－は電子の電気伝導度であり、ψ_ｅ ^－は電子の電位である。式７において、σ_Ｈ ^＋はプロトンの電気伝導度であり、ψ_Ｈ ^＋はプロトンの電位である。式６及び式７において、Ｒ_ｃａｔは、電子の触媒層５２における化学反応により生成された電流の密度であって、下記式８に示すバトラー・ボルマー方程式により得られる。

【0101】

【数8】

【0102】

式８において、ζ_ｃａｔは、触媒層５２における白金の有効面積であり、ｊ_ｃａｔは、参照交換電流密度である。［Ｏ_２］は酸素のモル分率であり、上記式５のガスの分布▽Ｘ_ｉである。［Ｏ_２］_ｒｅｆは参照の酸素モル分率である。γ_ｃａｔは、濃度次数であり、ａ_ｃａｔは、移動係数であり、Ｆはファラデー定数であり、Ｒは気体定数であり、Ｔは温度である。

【0103】

＜相変化を伴う現象のシミュレーション方法＞
相変化を伴う現象のシミュレーション方法は、例えば、シミュレーション装置１０によって図９に示すフローチャートに沿って実行される。図９に示すフローチャートでは、図７に示すステップＳ２０及びＳ６の処理に加えて、ステップＳ７－Ｓ１０の処理をさらに実行する。このステップＳ７の処理をステップＳ２０とＳ６との間に実行し、ステップ８－Ｓ１０の処理をＳ６の後に実行する。

【0104】

なお、図９に示すフローチャートでは、ステップＳ２０の処理の前に、図１１のステップＳ２１－Ｓ２３の処理、又は、図１２のステップＳ２４の処理を行ってもよい。また、図９のステップＳ２０の処理の前に、図１３のステップＳ２５の処理、又は、図１４のステップＳ２６の処理を行ってもよい。

【0105】

シミュレーション装置１０の演算処理部１１は、多孔質体２０の単位体積当たりの相変化速度（単位体積、単位時間当たりの相変化量ｒ_ｗ）、を算出する（ステップＳ５）。ここでは、ステップＳ２０の処理により液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗが求められる。

【0106】

また、演算処理部１１は、ガスの実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆを算出する（ステップＳ７）。この実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆは、例えば、上述の通り、ボルツマン方程式の考え方に基づくシミュレーション手法のガス拡散シミュレーションによる結果及び式５のＳｔｅｆａｎ－Ｍａｘｗｅｌｌ方程式を用いて得られる。

【0107】

続いて、演算処理部１１は、燃料電池５０における液水分布を算出する(ステップＳ６)。ここでは、上記式４により液水２４の充填率ｓを算出することで燃料電池５０における液水分布が求まる。

【0108】

また、演算処理部１１は、燃料電池５０における式４の実効的な浸透係数、式６、式７の実効的な電気電導度、実効的なプロトン伝導度及び式１０の実効的な熱伝導係数を式９の離散化されたボルツマン方程式に基づいて算出する。式９において、αは位相空間座標の指標であり、ｆ_ｉはｉ種粒子の分布関数であり、Ω_ｉｊ ^αはｉ種及びｊ種の衝突項である。式１０において、ｊ_Ｔは熱流束であり、σ_Ｔ ^ｅｆｆは実効的な熱伝導度であり、Ｔは温度である。

【0109】

【数9】

【0110】

【数10】

【0111】

続いて、演算処理部１１は、燃料電池５０のガスの分布を算出する(ステップＳ８)。ここでは、ステップＳ６で求めた実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆを用いた上記式５において、燃料電池５０の各条件を代入する。これにより、ガスのモル分率を算出することで、燃料電池５０におけるガスの分布が求まる。

【0112】

そして、演算処理部１１は、燃料電池５０の電位の分布を算出する(ステップＳ９)。ここでは、上記式６～８により電位を算出することで、燃料電池５０における電位の分布が求まる。

【0113】

続いて、演算処理部１１は、燃料電池５０の電位の分布が収束したか否かを判定する（ステップＳ１０）。ステップＳ９で算出した電位の変化量が所定値よりも大きければ、電位が収束していないとして（ステップＳ１０：ＮＯ）、ステップＳ６の処理に戻り、ステップＳ６及びＳ８～Ｓ１０の処理を繰り返す。一方、電位の変化量が所定値以下であれば、電位が収束したとして（ステップＳ１０：ＹＥＳ）、シミュレーションを終了する。

【0114】

このような相変化を伴う現象のシミュレーション方法では、液水２４の単位体積、単位時間当たりの相変化量ｒ_ｗに基づいて燃料電池５０における液水分布を算出する。また、細孔２２において液水２４以外の空隙２３に拡散するガスの実効的な相互拡散係数Ｄ_ｉｊ ^ｅｆｆを、ボルツマン方程式の考え方に基づくシミュレーション手法を用いて算出し、これに基づいて燃料電池５０におけるガスの分布を算出する。そして、燃料電池５０における液水分布及びガスの分布に基づいて、燃料電池５０の電位の分布を算出している。

【0115】

このように、液水分布を考慮することにより、燃料電池５０における過剰な液水２４によるフラッディング等の現象を踏まえて、燃料電池５０の電位を高精度にシミュレーションすることができる。

【0116】

また、上記ステップＳ６、Ｓ８及びＳ９において液水分布が求まり、この分布に対して定まるガスの分布におけるガスの反応によって、燃料電池５０の電位が求められる。このように、液水分布を考慮することにより、燃料電池５０における過剰な液水２４によるフラッディング等の現象を踏まえて、燃料電池５０の電位を高精度にシミュレーションすることができる。

【0117】

なお、上述した実施の形態で示した構成は一例であり、適宜変更してもよい。

【0118】

また、上記説明から、当業者にとっては、本発明の多くの改良や他の実施形態が明らかである。従って、上記説明は、例示としてのみ解釈されるべきであり、本発明を実行する最良の態様を当業者に教示する目的で提供されたものである。本発明の精神を逸脱することなく、その構造及び／又は機能の詳細を実質的に変更できる。

【産業上の利用可能性】

【0119】

本発明の相変化を伴う現象のシミュレーション方法は、高精度な相変化を伴う現象のシミュレーション方法等として有用である。

【符号の説明】

【0120】

１０ ：シミュレーション装置
１１ ：演算処理部
２０ ：多孔質体
２４ ：液水
５０ ：燃料電池

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】