(57)【要約】【課題】高精度な相変化を伴う現象のシミュレーション方法を提供する。【解決手段】相変化を伴う現象のシミュレーション方法では、多孔質体２０における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量を算出する。【選択図】図１０
2019216085.tif 000002
【特許請求の範囲】
【請求項１】
多孔質体における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量を算出する、相変化を伴う現象のシミュレーション方法。
【請求項２】
前記相変化量の算出前に、
前記多孔質体の構造情報を取得し、
前記多孔質体の構造情報に基づき、前記多孔質体における単位体積当たりの気液界面の面積を算出する、請求項１記載の相変化を伴う現象のシミュレーション方法。
【請求項３】
前記相変化量の算出前に、前記多孔質体における単位体積当たりの気液界面の面積を取得する、請求項１に記載の相変化を伴う現象のシミュレーション方法。
【請求項４】
前記相変化量の算出前に、分子動力学法を用いて前記気液界面の相変化速度を算出する、請求項１〜３のいずれか一項に記載の相変化を伴う現象のシミュレーション方法。
【請求項５】
前記相変化量の算出前に、前記気液界面の相変化速度を取得する、請求項１〜３のいずれか一項に記載の相変化を伴う現象のシミュレーション方法。
【請求項６】
前記多孔質体の単位体積、単位時間当たりの相変化量に基づいて、前記多孔質体を含む
燃料電池の液水分布を算出する、請求項１〜５のいずれか一項に記載の相変化を伴う現象のシミュレーション方法。
【請求項７】
前記多孔質体において液水以外の空隙に拡散するガスの実効的な相互拡散係数を、ボル
ツマン方程式の考え方に基づくシミュレーション手法を用いて算出し、
前記ガスの実効的な相互拡散係数に基づいて、前記燃料電池における前記ガスの分布を算出し、前記燃料電池における前記液水分布及び前記ガスの分布に基づいて、前記燃料電池における電位の分布を算出する、請求項６に記載の相変化を伴う現象のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、相変化を伴う現象のシミュレーション方法に関する。
【背景技術】
【０００２】
多孔質体を含む燃料電池では、多孔質体における液体の水（液水）の輸送により、その性能が影響を受ける。このため、例えば、特許文献１の燃料電池のシミュレーション方法では、触媒層の構造を算出し、水の輸送等に基づき触媒層の性能を算出し、触媒層の性能を基に燃料電池の発電性能を算出することにより、電解質膜と電極とを有する膜−電極を備える燃料電池の発電性能をシミュレーションしている。
【０００３】
また、多孔質体における液水の輸送は、液水に加わる圧力勾配に加え、蒸発及び凝縮の相変化にも影響を受ける。この液水の蒸発量及び凝縮量の相変化量について、例えば、非特許文献１の方法によれば、多孔質体における液水の充填率（飽和度）に比例するように決定されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１５−１３５７９０号公報
【非特許文献】
【０００５】
【非特許文献１】Ｅｎｅｒｇｙ    Ｖｏｌ．  ７３    ｐｐ６１８−６３４  （２０１４）
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明は、高精度な相変化を伴う現象のシミュレーション方法の構築を課題としている。
【課題を解決するための手段】
【０００７】
本発明のある態様に係る相変化を伴う現象のシミュレーション方法では、多孔質体における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量を算出する。
【発明の効果】
【０００８】
本発明は、高精度な相変化を伴う現象のシミュレーション方法を提供することができるという効果を奏する。
【０００９】
本発明の上記目的、他の目的、特徴、及び利点は、添付図面参照の下、以下の好適な実施態様の詳細な説明から明らかにされる。
【図面の簡単な説明】
【００１０】
【図１】本発明の実施の形態１に係る相変化を伴う現象のシミュレーション装置を示す機能ブロック図である。
【図２】相変化を伴う現象のシミュレーション方法に用いられる多孔質体を模式的に示した断面図である。
【図３】気液界面の相変化速度に関する解析領域を模式的に示す斜視図である。
【図４】気体分子の数の経時変化を示すグラフである。
【図５】その他の変形例に係るシミュレーション方法の一例を示すフローチャートである。
【図６】燃料電池を模式的に示す斜視図である。
【図７】本発明の実施の形態２に係る相変化を伴う現象のシミュレーション方法の一例を示すフローチャートである。
【図８】ガスの相互拡散係数に関する解析領域を模式的に示す斜視図である。
【図９】本発明の実施の形態３に係る相変化を伴う現象のシミュレーション方法の一例を示すフローチャートである。
【図１０】実施の形態１に係るシミュレーション方法の一例を示すフローチャートである。
【図１１】変形例１に係るシミュレーション方法の一例を示すフローチャートである。
【図１２】変形例２に係るシミュレーション方法の一例を示すフローチャートである。
【図１３】変形例３に係るシミュレーション方法の一例を示すフローチャートである。
【図１４】変形例４に係るシミュレーション方法の一例を示すフローチャートである。
【発明を実施するための形態】
【００１１】
（本発明の基礎となった知見）
本発明者らは、多孔質体における液水に関するシミュレーションの精度向上について検討した。液水の相変化は、液水の気液界面を介して水分子が移動することにより生じる。しかしながら、非特許文献１の方法では、液水の相変化量は、多孔質体における液水の充填率に基づいており、多孔質体における気液界面の面積が液水の相変化量に反映されていない。そこで、本発明者等は、多孔質体における気液界面の面積を考慮して液水の相変化量を求めることにより、シミュレーションの精度向上が図られることを見出した。
【００１２】
第１の発明に係る相変化を伴う現象のシミュレーション方法では、多孔質体における単位体積当たりの気液界面の面積、及び、前記気液界面の相変化速度から、前記多孔質体の単位体積、単位時間当たりの相変化量を算出する。
【００１３】
この方法によれば、液水の相変化量（蒸発量及び凝縮量）は、多孔質体における単位体積当たりの気液界面の面積に基づいて算出される。このため、液水の相変化量について、相変化が生じる気体と液水の気液界面を考慮した高精度なシミュレーションを行うことができる。
【００１４】
第２の発明に係る相変化を伴う現象のシミュレーション方法では、第１の発明において、前記相変化量の算出前に、前記多孔質体の構造情報を取得し、前記多孔質体の構造情報に基づき、前記多孔質体における単位体積当たりの気液界面の面積を算出する。この方法によれば、多孔質体の構造を考慮した液水の相変化量のシミュレーションを行うことができる。
【００１５】
第３の発明に係る相変化を伴う現象のシミュレーション方法では、第１の発明において、前記相変化量の算出前に、前記多孔質体における単位体積当たりの気液界面の面積を取得する。この方法によれば、取得した気液界面の面積を考慮した液水の相変化量を高精度にシミュレーションすることができる。
【００１６】
第４の発明に係る相変化を伴う現象のシミュレーション方法では、第１〜３のいずれかの発明において、前記相変化量の算出前に、分子動力学法を用いて前記気液界面の相変化速度を算出する。この方法によれば、単位体積当たりの単位面積当たりに蒸発する液水量を決定する相変化速度係数を決定することができる。この方法によれば、分子動力学法を用いて算出した相変化速度を考慮した液水の相変化量を高精度にシミュレーションすることができる。また、このシミュレーション方法においては、分子動力学法を用いることにより、実験等に比べて容易に気液界面の相変化速度を取得することができる。
【００１７】
第５の発明に係る相変化を伴う現象のシミュレーション方法では、第１〜３のいずれかの発明において、前記相変化量の算出前に、前記気液界面の相変化速度を取得する。この方法によれば、取得した相変化速度を考慮した液水の相変化量を高精度にシミュレーションすることができる。
【００１８】
第６の発明に係る相変化を伴う現象のシミュレーション方法では、第１〜５のいずれかの発明において、前記多孔質体の単位体積、単位時間当たりの相変化量に基づいて、前記多孔質体を含む燃料電池の液水分布を算出する。この方法によれば、燃料電池における液水分布について、液水の気液界面を考慮した高精度なシミュレーションを行うことができる。
【００１９】
第７の発明に係る相変化を伴う現象のシミュレーション方法では、第６の発明において、前記多孔質体において液水以外の空隙に拡散するガスの実効的な相互拡散係数を、ボルツマン方程式の考え方に基づくシミュレーション手法を用いて算出し、前記ガスの実効的な相互拡散係数に基づいて、前記燃料電池における前記ガスの分布を算出し、前記燃料電池における前記液水分布及び前記ガスの分布に基づいて、前記燃料電池における電位の分布を算出する。
【００２０】
この方法によれば、燃料電池の電位分布（発電性能予測）のシミュレーションを高精度に行うことができる。なお、ボルツマン方程式の考え方に基づくシミュレーション手法とは、格子ボルツマン法及び格子ガス法のことを意味している。
【００２１】
以下、本発明の実施形態を、図面を参照しながら具体的に説明する。なお、以下では全ての図面を通じて同一又は相当する要素には同一の参照符号を付して、その重複する説明を省略する。
【００２２】
（実施の形態１）
＜シミュレーション装置の構成＞
本発明の実施の形態１に係る相変化を伴う現象のシミュレーション装置１０は、図１に示すように、多孔質体における液体の水（液水）の相変化量等を予測する装置である。この多孔質体には、燃料電池セルに用いられるガス拡散層及び触媒層等が例示される。このため、シミュレーション装置１０は、例えば、燃料電池におけるシミュレーション装置として用いられるが、これ以外のシミュレーション装置にも用いることができる。
【００２３】
液水の相変化は、水の液体から気体への変化（蒸発）、及び、水の気体から液体への変化（凝縮）である。また、液水の相変化量は、蒸発する水の量（蒸発量）、又は、凝縮する水の量（凝縮量）である。
【００２４】
シミュレーション装置１０は、演算処理部１１及び記憶部１２を有しており、例えば、コンピュータによって構成されている。例えば、演算処理部１１としてコンピュータのＣＰＵ等が例示され、記憶部１２として演算処理部１１によりアクセス可能なメモリ等が例示される。
【００２５】
記憶部１２には、多孔質体における液水の相変化等をシミュレーションするプログラム及びこれに必要な情報が記憶されている。なお、このプログラム及びこれに必要な情報は、コンピュータに内蔵されたメモリに記憶されたものに限定されず、他の記憶媒体に記憶されたものであってもよいし、入力装置によって入力されたものであってもよいし、ネットワークを介して受信されたものであってもよい。
【００２６】
演算処理部１１は、記憶部１２に記憶されたプログラム等のソフトウェアを読み出して実行することにより、多孔質体における液水の相変化等をシミュレーションする。シミュレーション装置１０は、それぞれ、単独の装置で構成されていても構わないし、互いに協働する複数の装置で構成されていても構わない。
【００２７】
＜多孔質体の構成＞
シミュレーションの対象とする図２の多孔質体２０について説明する。多孔質体２０には、壁部分２１及び多数の細孔２２が設けられている。壁部分２１は、例えば、樹脂及びカーボン担体等の有機物、ガラス等の無機物、並びに、その混合物等により形成されている。
【００２８】
細孔２２は、多孔質体２０において周囲が壁部分２１により取り囲まれていた空間である。細孔２２に液体状の水（液水２４）が分布している場合、細孔２２から液水２４を除いた空間（空隙２３）が、ガスの粒子が移動可能な空間となる。このため、空隙２３は、多孔質体２０において壁部分２１及び液水２４以外の領域である。例えば、燃料電池の場合、ガスの粒子としては、水素、酸素及び窒素等が例示される。
【００２９】
＜相変化を伴う現象のシミュレーション方法＞
相変化を伴う現象のシミュレーション方法は、例えば、図１のシミュレーション装置１０によって図１０に示すフローチャートに沿って実行される。ここでは、図２に示す多孔質体２０の単位体積及び単位時間当たりの液水２４の相変化量についてシミュレーションする。
【００３０】
シミュレーション装置１０の演算処理部１１は、多孔質体２０における単位体積当たりの気液界面の面積、及び、気液界面の相変化速度から、多孔質体２０の単位体積、単位時間当たりの相変化量を算出する（ステップＳ２０）。
【００３１】
この液水２４の相変化は、多孔質体２０の細孔２２において空隙２３に接する液水２４の気液界面（液水２４の表面）において生じる。多孔質体２０における単位体積当たりの気液界面の面積（比表面積Ｓｆ）は、単位体積の多孔質体２０において細孔２２の空隙２３に接している液水２４の気液界面の面積である。また、気液界面の相変化速度（相変化速度係数ｋ）は、気液界面において液水２４が蒸発する速度、又は、気液界面において液水２４に凝縮する速度である。
【００３２】
続いて、演算処理部１１は、相変化速度係数ｋ及び比表面積Ｓｆに基づき、下記式３から多孔質体２０の単位体積及び単位時間当たりの液水２４の相変化量（多孔質体２０の単位体積当たりの気液界面の相変化速度（単位体積、単位時間当たりの相変化量ｒｗ））を算出する。
【００３３】
この式３において、ｋは、相変化速度係数である。液水２４の水蒸気分圧ｐｗｖ、液水２４の飽和蒸気圧ｐｓａｔ、温度Ｔは、求めたいシミュレーション条件により定められる。
【００３４】
【数１】
2019216085.tif 000003
【００３５】
このような相変化を伴う現象のシミュレーション方法では、液水２４の比表面積Ｓｆ及び相変化速度係数ｋから、単位体積、単位時間当たりの相変化量ｒｗを算出している。このように、液水２４の比表面積Ｓｆによって、多孔質体２０において相変化が生じる液水２４の気液界面の面積（表面積）を考慮して、液水２４の単位体積、単位時間当たりの相変化量ｒｗを高精度にシミュレーションすることができる。
【００３６】
＜変形例１＞
変形例１に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、多孔質体２０の構造情報を取得し、多孔質体２０の構造情報に基づき、多孔質体２０における単位体積当たりの気液界面の面積を算出する。このシミュレーション方法は、図１のシミュレーション装置１０によって図１１に示すフローチャートに沿って実行される。
【００３７】
具体的には、演算処理部１１は、図１１のシミュレーション方法において対象とする多孔質体２０の構造情報を取得する（ステップＳ２１）。この多孔質体２０の構造情報は、多孔質体２０の形状及び材質に関する情報を含む。
【００３８】
例えば、多孔質体２０の壁部分２１の材質情報としては、壁部分２１に対する液水２４の接触角θが例示され、これは実験等より予め定めされうる。また、多孔質体２０内の液水２４の情報として、表面張力γがあり、これは実験等より予め定めされている。
【００３９】
壁部分２１の形状情報は、多孔質体２０の画像情報等により得られる。例えば、ＦＩＢ−ＳＥＭ（Ｆｏｃｕｓｅｄ  Ｉｏｎ  Ｂｅａｍ−Ｓｃａｎｎｉｎｇ  Ｅｌｅｃｔｒｏｎ  Ｍｉｃｒｏｓｃｏｐｅ）を用いて、多孔質体２０について所定の方向において異なる位置で複数の断面画像を連続的に撮影する。
【００４０】
この複数の断面画像を積層した画像（積層画像）を画像処理によって二値化の作業を行う。これにより、細孔面２０ａの位置情報を含むシミュレーション用の計算メッシュが、多孔質体２０の形状情報として作成される。ここで、積層画像に基づいた多孔質体２０の形状情報の取得には、Ｍａｔｈ２Ｍａｒｋｅｔ社製のＧｅｏＤｉｃｔ（ドイツの登録商標）等の解析ソフトウェアを用いてもよい。
【００４１】
そして、演算処理部１１は、ステップＳ２１において情報が取得された多孔質体２０における液水２４の分布（液水分布）を取得する（ステップＳ２２）。液水分布は、多孔質体２０の細孔２２に配置される液水２４の位置情報（座標等）であって、以下のとおり、式２及びポアモルフォロジー法による液水２４を詰め込むシミュレーションを行うことにより求められる。
【００４２】
【数２】
2019216085.tif 000004
【００４３】
例えば、これは、Ｍａｔｈ２Ｍａｒｋｅｔ社製の解析ソフトＧｅｏＤｉｃｔ等のソフトウェアに用いることにより算出することが可能である。具体的には、まず、多孔質体２０の壁部分２１の材質情報、及び、圧力（毛管圧）により、多孔質体２０内に液水２４が存在可能な空隙の半径をヤング・ラプラスの式（式２）により決定する。式２において、ｐｃは毛管圧［Ｐａ］である。γは液水２４の表面張力［Ｎ／ｍ］であり、θは壁部分２１に対する液水２４の接触角［°］であって、これらは壁部分２１の材質情報から求められる。ｒは、空隙の半径［ｍ］である。
【００４４】
多孔質体２０における壁部分２１に対する液水２４の接触角θが９０°を超える（θ＞９０°）場合、壁部分２１は撥水性である。このため、液水２４は、毛管圧Ｐｃによって細孔２２から押し出される力を受ける。毛管圧Ｐｃから式２に基づいて液水２４が浸入できる空隙２３の最小の半径ｒを求める。そして、多孔質体２０の所定面から液水２４を順次、詰め込むシミュレーションを、空隙２３の最小の半径ｒを基にポアモルフォロジー法により実行する。これにより、多孔質体２０に配置される液水２４の位置情報が液水分布として求められる。
【００４５】
一方、液水２４の接触角θが９０°を下回る（θ＜９０°）場合、多孔質体２０の壁部分２１は親水性である。このため、液水２４は、毛管圧Ｐｃによって細孔２２に取り込まれる力を受ける。毛管圧Ｐｃから式２に基づいて気体により液水２４が排出される空隙２３の最大の半径ｒを求める。そして、上記多孔質体２０の形状情報に基づく多孔質体２０の所定面から液水２４が記載により順次、押し出されるシミュレーションを、空隙２３の最大の半径ｒを基に上記ポアモルフォロジー法により実行する。これにより、多孔質体２０に配置される液水２４の位置情報が液水分布として求められる。
【００４６】
ここで、毛管圧Ｐｃを変化させることにより、多孔質体２０における液水２４の充填率が変化する。このため、式２において毛管圧Ｐｃを変えながら、シミュレーションを実行することにより、液水分布が充填率の関数として得られる。なお、液水２４の充填率は０から１の範囲に設定される。液水２４の充填率が０は、細孔２２に液水２４が配置されていない状態を示し、液水２４の充填率が１は、細孔２２に液水２４が満たされている状態を示す。
【００４７】
例えば、Ｍａｔｈ２Ｍａｒｋｅｔ社製の解析ソフトＧｅｏＤｉｃｔを用いることにより、多孔質体２０の構造情報及び液水２４の充填率から液水分布が決定される。
【００４８】
続いて、演算処理部１１は、ステップＳ２２の液水分布に基づき、多孔質体２０の単位体積当たりの気液界面の面積（比表面積Ｓｆ）を算出する（ステップＳ２３）。
【００４９】
ここで、Ｓ２１における壁部分２１の形状情報から、細孔２２の位置が求められる。また、Ｓ２２における液水分布から、細孔２２における液水２４の位置が求められる。このため、これらの位置から、多孔質体２０における空隙２３に接する液水２４の気液界面の面積（表面積）を算出する。
【００５０】
そして、気液界面の面積を多孔質体２０の体積により割ることにより、液水２４の比表面積Ｓｆが算出される。ここで、液水分布が充填率の関数で表されている場合、これに基づく液水２４の比表面積Ｓｆも充填率の関数で定められる。例えば、液水２４の比表面積Ｓｆは、Ｍａｔｈ２Ｍａｒｋｅｔ社製のＧｅｏＤｉｃｔ等のソフトウェアを用いることにより得られる。
【００５１】
続いて、演算処理部１１は、気液界面の相変化速度係数ｋ、及び、ステップＳ２３の液水２４の比表面積Ｓｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒｗを算出する（ステップＳ２０）。
【００５２】
このように、多孔質体２０の構造情報に基づいて算出した液水２４の比表面積Ｓｆを用いて、多孔質体２０の単位体積、単位時間当たりの相変化量ｒｗを算出している。これにより、多孔質体２０の構造を考慮した、液水２４の単位体積、単位時間当たりの相変化量ｒｗを高精度にシミュレーションすることができる。
【００５３】
＜変形例２＞
変形例２に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、多孔質体２０における単位体積当たりの気液界面の面積を取得する。このシミュレーション方法は、図１１に示すフローチャートにおけるステップＳ２１〜Ｓ２３の処理に代えて、図１２に示すフローチャートにおけるステップＳ２４の処理が実行される。
【００５４】
具体的には、演算処理部１１は、図１２のシミュレーション方法において対象とする多孔質体２０の単位体積当たりの気液界面の面積（比表面積Ｓｆ）を取得する（ステップＳ２４）。この比表面積Ｓｆは、シミュレーション装置１０に対する入力値であってもよく、例えば、記憶部１２に記憶されたものに限定されず、他の記憶媒体に記憶されたものであってもよいし、入力装置によってシミュレーション装置１０に入力されたものであってもよいし、ネットワークを介して受信されたものであってもよい。
【００５５】
そして、演算処理部１１は、気液界面の相変化速度係数ｋ、及び、ステップＳ２４の液水２４の比表面積Ｓｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒｗを算出する（ステップＳ２０）。
【００５６】
このように、取得した液水２４の比表面積Ｓｆを用いて、多孔質体２０の単位体積、単位時間当たりの相変化量ｒｗを算出する。これにより、取得した比表面積Ｓｆに応じた、液水２４の単位体積、単位時間当たりの相変化量ｒｗを高精度にシミュレーションすることができる。
【００５７】
＜変形例３＞
変形例３に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、分子動力学法を用いて気液界面の相変化速度を算出する。このシミュレーション方法は、図１のシミュレーション装置１０によって図１３に示すフローチャートに沿って実行される。
【００５８】
具体的には、演算処理部１１は、気液界面の相変化速度（相変化速度係数ｋ）を算出する（ステップＳ２５）。この相変化の速度（相変化速度）は、例えば、図３に示す解析領域３０において分子動力学法によるシミュレーションを行うことによって求められる。
【００５９】
解析領域３０は、例えば、直方体形状の空間である。この解析領域３０に、膜のように薄い直方体形状の液水２４（液膜３１）が、配置されている。液膜３１は、解析領域３０の長手方向に対して直交する方向に拡がり、解析領域３０を長手方向に２つの部分（第１空間３２、第２空間３３）に区分けするように配置されている。このため、液膜３１は、第１空間３２に接する一方面（第１主面３１ａ）、及び、第２空間３３に接する他方面（第２主面３１ｂ）を有している。第１主面３１ａ及び第２主面３１ｂから成る液水２４の気液界面において液膜３１の水分子が相変化する。
【００６０】
解析領域３０において平衡状態の液膜３１に熱エネルギーを加える。これにより、液膜３１の水分子が液水２４の気液界面から蒸発し、気体分子となる。この気体分子の数をカウントする。この気体分子の数の経時変化が図４のグラフに求められる。
【００６１】
図４のグラフに示すように、時刻０から時刻ａまでは気体分子の数の経時変化が所定値以下であり、液膜３１が平衡状態である。そして、時刻ａで所定量のエネルギーを液膜３１に加える。この時刻ａにおいて気体分子の数の経時変化が増加し、水分子の蒸発が開始する。そして、時刻ｂにおいて気体分子の数の経時変化が減少し、水分子の蒸発が終了して、液膜３１が平衡状態となる。この時刻ａから時刻ｂまでの気体分子の数の経時変化から、気液界面の相変化速度ｒ０が求められる。
【００６２】
この気液界面の相変化速度ｒ０及びシミュレーション条件を下記式１に代入することにより、気液界面の単位面積当たりの液水２４の相変化速度（相変化速度係数ｋ）が算出される。この式１において、Ｓ０は、単位体積当たりの液水２４の気液界面の面積（液水２４の表面積）であって、液膜３１の第１主面３１ａ及び第２主面３１ｂの合計面積を解析領域３０の体積で割った値である。ｐｗｖは、液水２４の水蒸気分圧であり、ｐｓａｔは、液水２４の飽和蒸気圧であり、Ｔは、温度であって、これらはシミュレーション条件により定められる。Ｒは、気体定数であり、ＭＨ２Ｏは、液水２４の分子量である。
【００６３】
【数３】
2019216085.tif 000005
【００６４】
そして、演算処理部１１は、ステップＳ２５の気液界面の相変化速度係数ｋ、及び、液水２４の比表面積Ｓｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒｗを算出する（ステップＳ２０）。
【００６５】
このように、算出した気液界面の相変化速度を用いて、多孔質体２０の単位体積、単位時間当たりの相変化量ｒｗを算出している。これにより、算出条件に応じた、液水２４の単位体積、単位時間当たりの相変化量ｒｗを高精度にシミュレーションすることができる。
【００６６】
また、このシミュレーション方法においては、分子動力学法によるシミュレーションによって気液界面の相変化速度を取得可能である。この相変化速度から気液界面の相変化速度係数ｋを実験等に比べて容易に取得することができる。
【００６７】
なお、ステップＳ２０において、液水２４の比表面積Ｓｆは、例えば、変形例１のようにシミュレーション装置１０によって算出されてもよいし、変形例２のようにシミュレーション装置１０によって取得されてもよい。
【００６８】
例えば、気液界面の相変化速度係数ｋ、及び、液水２４の比表面積Ｓｆがシミュレーション装置１０によって算出される場合、相変化を伴う現象のシミュレーション方法は、図５に示すフローチャートに沿って実行される。この図５のステップＳ１は図１３のステップＳ２５に相当し、図５のステップＳ２〜Ｓ４は図１１のステップＳ２１〜Ｓ２３にそれぞれ相当し、図５のステップＳ５は図１１のステップＳ２０に相当する。
【００６９】
＜変形例４＞
変形例４に係る相変化を伴う現象のシミュレーション方法では、相変化量の算出前に、気液界面の相変化速度を取得する。このシミュレーション方法は、図１３に示すフローチャートにおけるステップＳ２５の処理に代えて、図１４に示すフローチャートにおけるステップＳ２６の処理が実行される。
【００７０】
具体的には、演算処理部１１は、気液界面の相変化速度（相変化速度係数ｋ）を取得する（ステップＳ２５）。この相変化速度係数ｋは、シミュレーション装置１０に対する入力値であってもよく、例えば、記憶部１２に記憶されたものに限定されず、他の記憶媒体に記憶されたものであってもよいし、入力装置によって入力されたものであってもよいし、ネットワークを介して受信されたものであってもよい。
【００７１】
そして、演算処理部１１は、ステップＳ２６の気液界面の相変化速度係数ｋ、及び、液水２４の比表面積Ｓｆに基づき、上記式３から多孔質体２０の単位体積、単位時間当たりの相変化量ｒｗを算出する（ステップＳ２０）。
【００７２】
このように、取得した相変化速度係数ｋに応じた、液水２４の単位体積、単位時間当たりの相変化量ｒｗを高精度にシミュレーションすることができる。
【００７３】
なお、ステップＳ２０において、液水２４の比表面積Ｓｆは、例えば、変形例１のようにシミュレーション装置１０によって算出されてもよいし、変形例２のようにシミュレーション装置１０によって取得されてもよい。
【００７４】
（実施の形態２）
本発明の実施の形態２に係るシミュレーション装置１０は、図６に示す燃料電池５０における液水２４の分布（液水分布）をその発電性能として予測する装置である。
【００７５】
燃料電池５０は、例えば、電解質膜５１、これを互いの間に挟む一対の触媒層５２、これを互いの間に挟む一対のガス拡散層５３、及び、これを互いの間に挟むガス流路５４により構成されている。このうち、触媒層５２及びガス拡散層５３の少なくともいずれか１つの部材が多孔質体２０により形成されている。
【００７６】
シミュレーション装置１０による相変化を伴う現象のシミュレーション方法では、燃料電池５０内における液水２４の充填率を算出する。
【００７７】
＜燃料電池における液水分布の決定＞
多孔質体２０の単位体積当たりの相変化速度（単位体積、単位時間当たりの相変化量ｒｗ）に基づいて、燃料電池５０における液水分布を下記式４により算出することができる。ここで、燃料電池５０における条件を式４に代入する。
【００７８】
式４において、εは細孔２２に対する空隙２３の割合（空隙率）であり、ρｌは液水２４の密度であり、ｓは細孔２２内の液水２４の充填率であり、ｔは時間［ｓｅｃ］である。また、Ｋは液水２４の透過率であり、ｒは液水２４の充填率の次数である。Ｐｃは毛管圧であり、上記式１により求められる。μｌは液水２４の粘性係数である。ｒｗは液水２４の単位体積、単位時間当たりの相変化量であり、上記式３から求められる。
【００７９】
【数４】
2019216085.tif 000006
【００８０】
このように、燃料電池５０内での液水２４の充填率ｓの分布が式４により求められる。これは、例えば、ＦＬＵＥＮＴ（登録商標）等のソフトウェアに導入することにより算出される。
【００８１】
＜相変化を伴う現象のシミュレーション方法＞
相変化を伴う現象のシミュレーション方法は、例えば、シミュレーション装置１０によって図７に示すフローチャートに沿って実行される。図７に示すフローチャートでは、図１０に示すステップＳ２０の処理に加えて、ステップＳ６の処理をさらに実行する。なお、図７のステップＳ２０の処理の前に、図１１のステップＳ２１〜Ｓ２３の処理、又は、図１２のステップＳ２４の処理を行ってもよい。また、図７のステップＳ２０の処理の前に、図１３のステップＳ２５の処理、又は、図１４のステップＳ２６の処理を行ってもよい。
【００８２】
シミュレーション装置１０の演算処理部１１は、多孔質体２０の単位体積当たりの相変化速度（単位体積、単位時間当たりの相変化量ｒｗ）を算出する（ステップＳ２０）。ここでは、ステップＳ２０の処理により液水２４の単位体積、単位時間当たりの相変化量ｒｗが求められる。
【００８３】
また、演算処理部１１は、液水２４の単位体積、単位時間当たりの相変化量ｒｗに基づいて、燃料電池５０における液水２４の分布（液水分布）を算出する（ステップＳ６）。この燃料電池５０における液水分布は、上述の通り、式４により求められる。
【００８４】
この相変化を伴う現象のシミュレーション方法では、液水２４の単位体積、単位時間当たりの相変化量ｒｗに基づいて燃料電池５０における液水分布を算出している。これにより、液水２４の比表面積Ｓｆに基づいた液水２４の単位体積、単位時間当たりの相変化量ｒｗによって、多孔質体２０において相変化が生じる液水２４の気液界面の面積を考慮して、燃料電池５０における液水分布を高精度にシミュレーションすることができる。
【００８５】
また、式４における液水２４の単位体積、単位時間当たりの相変化量ｒｗに、多孔質体２０における液水２４の輸送が依存する。このため、多孔質体２０から成る部材を含む燃料電池５０において、液水２４の輸送を考慮した高精度のシミュレーションを行うことができる。
【００８６】
（実施の形態３）
本発明の実施の形態３に係るシミュレーション装置１０は、図６に示す燃料電池５０の電位の分布をその発電性能として予測する装置である。燃料電池５０の電位分布は、燃料電池５０における液水２４の分布（液水分布）及びガスの分布に基づいて求められる。燃料電池５０におけるガスの分布は、ガスの実効的な相互拡散係数に基づいて算出される。
【００８７】
＜ガスの実効的な相互拡散係数の決定＞
ガスの粒子が細孔２２における空隙２３において移動しながら、ガスが拡散する。このガスの拡散について、ボルツマン方程式の考え方に基づくシミュレーション手法を用いたガス拡散シミュレーションでは、図８に示す流路構造の解析領域４０が用いられる。
【００８８】
解析領域４０は、多孔質体２０、第１流路４１及び第２流路４２を有している。第１流路４１及び第２流路４２は、例えば、断面が矩形状であって、第１方向に長く延びており、互いに間隔を空けて互いに平行に設けられている。多孔質体２０は、例えば、直方体形状であって、第１流路４１と第２流路４２との間に配置されており、第１方向に直交する第２方向に延びている。
【００８９】
第１方向における多孔質体２０の寸法は第１流路４１及び第２流路４２の各寸法よりも小さい。第２方向における多孔質体２０の寸法は第１流路４１及び第２流路４２の各寸法よりも大きい。第１方向及び第２方向に直交する第３方向における孔質体の寸法は第１流路４１及び第２流路４２の各寸法と等しい。
【００９０】
このような解析領域４０において、多孔質体２０は第１流路４１及び第２流路４２に接続されており、多孔質体２０の細孔２２は第１流路４１及び第２流路４２とそれぞれ連通している。細孔２２は第１流路４１側の面と第２流路４２側の面との間に長く延びており、第１流路４１と第２流路４２とは細孔２２を介して互いに連通している。
【００９１】
ここで、液水２４が多孔質体２０に充填され、細孔２２に形成された空隙２３と第１流路４１及び第２流路４２とが連通する。この第１流路４１及び第２流路４２には互いに異なる種類のガス（第１ガス、第２ガス）を流入する。この流入するガスの流速は、実際の燃料電池の触媒層の中を流れるガスの平均的な流速に設定する。
【００９２】
第１ガスの粒子（第１ガス粒子ｉ）は第１流路４１の入口から流入し、その一部が多孔質体２０の空隙２３に流れ込んで拡散し、残りは第１流路４１の出口から流出する。第２ガスの粒子（第２ガス粒子ｊ）は第２流路４２の入口から流入し、その一部が多孔質体２０の空隙２３に流れ込んで拡散し、残りは第２流路４２の出口から流出する。
【００９３】
これにより、第１流路４１の出口における第１ガスのモル分率、及び、第２流路４２の出口における第２ガスのモル分率は、空隙２３における第１ガス及び第２ガスの拡散状態（ガスの相互拡散）に依存する。このため、上記格子ボルツマン法のガス拡散シミュレーションによる各流路の出口における第１ガスのモル分率及び第２ガスのモル分率が、下記式５で示されるＳｔｅｆａｎ−Ｍａｘｗｅｌｌ方程式のモル分率に一致するように、実効的な相互拡散係数Ｄｉｊｅｆｆを調整することによって、実効的な相互拡散係数Ｄｉｊｅｆｆが求められる。例えば、これは、ＦＬＵＥＮＴ等の熱流体解析ソフトウェアを利用することにより算出することが可能である。
【００９４】
【数５】
2019216085.tif 000007
【００９５】
式５において、Ｘｉは第１ガス粒子ｉのモル分率であり、Ｘｉは第２ガス粒子ｊのモル分率である。Ｖｉは第１ガス粒子ｉの速度であり、Ｖｊは第２ガス粒子ｊの速度である。Ｄｉｊｅｆｆは、第１ガス粒子ｉと第２ガス粒子ｊとの実効的な相互拡散係数（ガスの実効的な相互拡散係数）である。ガスの実効的な相互拡散係数Ｄｉｊｅｆｆは、多孔質体２０の形状及び多孔質体２０における液水分布を考慮したガスの相互拡散係数である。
【００９６】
この多孔質体２０における液水分布は、多孔質体２０における液水２４の充填率に依存する。よって、ガスの実効的な相互拡散係数Ｄｉｊｅｆｆは、充填率を変数とする関数として得られる。
【００９７】
＜燃料電池の電位分布の決定＞
燃料電池５０の電位は、下記式６、式７及び式８で表される電子及びプロトンに関する方程式により算出される。また、例えば、ＦＬＵＥＮＴ等のソフトウェアに用いることにより算出されうる。
【００９８】
【数６】
2019216085.tif 000008
【００９９】
【数７】
2019216085.tif 000009
【０１００】
式６において、σｅ−は電子の電気伝導度であり、ψｅ−は電子の電位である。式７において、σＨ＋はプロトンの電気伝導度であり、ψＨ＋はプロトンの電位である。式６及び式７において、Ｒｃａｔは、電子の触媒層５２における化学反応により生成された電流の密度であって、下記式８に示すバトラー・ボルマー方程式により得られる。
【０１０１】
【数８】
2019216085.tif 000010
【０１０２】
式８において、ζｃａｔは、触媒層５２における白金の有効面積であり、ｊｃａｔは、参照交換電流密度である。［Ｏ２］は酸素のモル分率であり、上記式５のガスの分布▽Ｘｉである。［Ｏ２］ｒｅｆは参照の酸素モル分率である。γｃａｔは、濃度次数であり、ａｃａｔは、移動係数であり、Ｆはファラデー定数であり、Ｒは気体定数であり、Ｔは温度である。
【０１０３】
＜相変化を伴う現象のシミュレーション方法＞
相変化を伴う現象のシミュレーション方法は、例えば、シミュレーション装置１０によって図９に示すフローチャートに沿って実行される。図９に示すフローチャートでは、図７に示すステップＳ２０及びＳ６の処理に加えて、ステップＳ７−Ｓ１０の処理をさらに実行する。このステップＳ７の処理をステップＳ２０とＳ６との間に実行し、ステップ８−Ｓ１０の処理をＳ６の後に実行する。
【０１０４】
なお、図９に示すフローチャートでは、ステップＳ２０の処理の前に、図１１のステップＳ２１−Ｓ２３の処理、又は、図１２のステップＳ２４の処理を行ってもよい。また、図９のステップＳ２０の処理の前に、図１３のステップＳ２５の処理、又は、図１４のステップＳ２６の処理を行ってもよい。
【０１０５】
シミュレーション装置１０の演算処理部１１は、多孔質体２０の単位体積当たりの相変化速度（単位体積、単位時間当たりの相変化量ｒｗ）、を算出する（ステップＳ５）。ここでは、ステップＳ２０の処理により液水２４の単位体積、単位時間当たりの相変化量ｒｗが求められる。
【０１０６】
また、演算処理部１１は、ガスの実効的な相互拡散係数Ｄｉｊｅｆｆを算出する（ステップＳ７）。この実効的な相互拡散係数Ｄｉｊｅｆｆは、例えば、上述の通り、ボルツマン方程式の考え方に基づくシミュレーション手法のガス拡散シミュレーションによる結果及び式５のＳｔｅｆａｎ−Ｍａｘｗｅｌｌ方程式を用いて得られる。
【０１０７】
続いて、演算処理部１１は、燃料電池５０における液水分布を算出する(ステップＳ６)。ここでは、上記式４により液水２４の充填率ｓを算出することで燃料電池５０における液水分布が求まる。
【０１０８】
また、演算処理部１１は、燃料電池５０における式４の実効的な浸透係数、式６、式７の実効的な電気電導度、実効的なプロトン伝導度及び式１０の実効的な熱伝導係数を式９の離散化されたボルツマン方程式に基づいて算出する。式９において、αは位相空間座標の指標であり、ｆｉはｉ種粒子の分布関数であり、Ωｉｊαはｉ種及びｊ種の衝突項である。式１０において、ｊＴは熱流束であり、σＴｅｆｆは実効的な熱伝導度であり、Ｔは温度である。
【０１０９】
【数９】
2019216085.tif 000011
【０１１０】
【数１０】
2019216085.tif 000012
【０１１１】
続いて、演算処理部１１は、燃料電池５０のガスの分布を算出する(ステップＳ８)。ここでは、ステップＳ６で求めた実効的な相互拡散係数Ｄｉｊｅｆｆを用いた上記式５において、燃料電池５０の各条件を代入する。これにより、ガスのモル分率を算出することで、燃料電池５０におけるガスの分布が求まる。
【０１１２】
そして、演算処理部１１は、燃料電池５０の電位の分布を算出する(ステップＳ９)。ここでは、上記式６〜８により電位を算出することで、燃料電池５０における電位の分布が求まる。
【０１１３】
続いて、演算処理部１１は、燃料電池５０の電位の分布が収束したか否かを判定する（ステップＳ１０）。ステップＳ９で算出した電位の変化量が所定値よりも大きければ、電位が収束していないとして（ステップＳ１０：ＮＯ）、ステップＳ６の処理に戻り、ステップＳ６及びＳ８〜Ｓ１０の処理を繰り返す。一方、電位の変化量が所定値以下であれば、電位が収束したとして（ステップＳ１０：ＹＥＳ）、シミュレーションを終了する。
【０１１４】
このような相変化を伴う現象のシミュレーション方法では、液水２４の単位体積、単位時間当たりの相変化量ｒｗに基づいて燃料電池５０における液水分布を算出する。また、細孔２２において液水２４以外の空隙２３に拡散するガスの実効的な相互拡散係数Ｄｉｊｅｆｆを、ボルツマン方程式の考え方に基づくシミュレーション手法を用いて算出し、これに基づいて燃料電池５０におけるガスの分布を算出する。そして、燃料電池５０における液水分布及びガスの分布に基づいて、燃料電池５０の電位の分布を算出している。
【０１１５】
このように、液水分布を考慮することにより、燃料電池５０における過剰な液水２４によるフラッディング等の現象を踏まえて、燃料電池５０の電位を高精度にシミュレーションすることができる。
【０１１６】
また、上記ステップＳ６、Ｓ８及びＳ９において液水分布が求まり、この分布に対して定まるガスの分布におけるガスの反応によって、燃料電池５０の電位が求められる。このように、液水分布を考慮することにより、燃料電池５０における過剰な液水２４によるフラッディング等の現象を踏まえて、燃料電池５０の電位を高精度にシミュレーションすることができる。
【０１１７】
なお、上述した実施の形態で示した構成は一例であり、適宜変更してもよい。
【０１１８】
また、上記説明から、当業者にとっては、本発明の多くの改良や他の実施形態が明らかである。従って、上記説明は、例示としてのみ解釈されるべきであり、本発明を実行する最良の態様を当業者に教示する目的で提供されたものである。本発明の精神を逸脱することなく、その構造及び／又は機能の詳細を実質的に変更できる。
【産業上の利用可能性】
【０１１９】
本発明の相変化を伴う現象のシミュレーション方法は、高精度な相変化を伴う現象のシミュレーション方法等として有用である。
【符号の説明】
【０１２０】
１０      ：シミュレーション装置
１１      ：演算処理部
２０      ：多孔質体
２４      ：液水
５０      ：燃料電池
【図１】
2019216085.tif 000013
【図２】
2019216085.tif 000014
【図３】
2019216085.tif 000015
【図４】
2019216085.tif 000016
【図５】
2019216085.tif 000017
【図６】
2019216085.tif 000018
【図７】
2019216085.tif 000019
【図８】
2019216085.tif 000020
【図９】
2019216085.tif 000021
【図１０】
2019216085.tif 000022
【図１１】
2019216085.tif 000023
【図１２】
2019216085.tif 000024
【図１３】
2019216085.tif 000025
【図１４】
2019216085.tif 000026