JP2010211356A - 流動場解析プログラムおよび電場解析プログラム - Google Patents

Abstract

【課題】 解析速度の向上や解析負荷の軽減をさせた流動場解析プログラムおよび電場解析プログラムを提供する。
【解決手段】 流動データ記憶部と、圧力項を除いた流体力学における運動量保存の式を時間積分して得られる式（１）を空間的に離散化し、流動場データ記憶部から読み出した流速データを代入して暫定的な流速を計算する暫定流速計算部と、運動量保存の式と質量保存の式とを連成し、かつ暫定的な流速に係る係数部にダルシー則抵抗を加えて得られる式（２）を空間的に離散化し、流動場データ記憶部から読み出した圧力データおよび前記暫定的な流速を代入して、圧力が収束するまで繰り返し計算をする圧力計算部と、式（３）に、前記暫定的な流速と圧力計算部により計算された圧力の値とを代入して所定時刻後の流速に補正する流速補正部としてコンピュータを機能させる。
【選択図】 図１

Description

本発明は、燃料電池等に利用可能な流動場解析プログラムおよび電場解析プログラムの解析速度の向上や解析負荷の軽減に関するものである。

燃料電池は、電力発生に伴うＣＯ_２排出のないクリーンな次世代エネルギーとして注目されている。この燃料電池は、多孔質材や電解質膜等から構成されており、水素等の燃料剤と酸素等の酸化剤を化学反応させることにより電力を取り出すシステムである。従来、この燃料電池の発電効率向上を目的として、燃料電池内の化学反応解析や流動解析等を行うプログラムの開発がなされている。

例えば、特開２００５−４４６０２号公報においては、燃料電池の物理的及び化学的特性の解析方法並びにプログラムが提案されている（特許文献１）。この特許文献１に記載のプログラムは、はじめに初期値の設定を行い、その後、酸素濃度、電流密度、水素濃度等、温度分布、化学反応をそれぞれ解析し、これらのデータをもとに燃料電池全体の電位が釣り合うか否かの判別をし、さらに電流値が実験値と一致するか否かの判別を行うものであり、電気化学的拘束条件下の熱と流れを容易に、かつ精度よく解析することができると記載されている。

特開２００５−４４６０２号公報

しかしながら、特許文献１に記載された発明を含む、これまでの燃料電池の数値解析プログラムは以下の問題点を有している。
１．燃料剤等の流動の解析において、多孔質材の内部では一般的に用いられる圧力や流速の時間変化が主に対流項に支配されているとする仮定ができないため、圧力計算の収束が悪くなる。
２．電気化学反応の解析において、解析結果が設定する初期値に依存してしまう。そのため、初期値によっては収束性が悪くなったり、解が得られない場合がある。
３．電気化学反応の解析において、解析精度が多孔質材と電解質膜との間にある薄い触媒層近傍の格子解像度に依存する。そのため、格子を十分に小さくする必要があり、計算負荷が大きくなる。

本発明は、このような問題点を解決するためになされたものであって、解析速度の向上や解析負荷の軽減をさせた流動場解析プログラムおよび電場解析プログラムを提供することを目的としている。

本発明に係る流動場解析プログラムは、流体力学における運動量保存の式と流体力学における質量保存の式とを連成させて所定の流動場の流速および圧力を計算する流動場解析プログラムであって、
流速データ（v）および圧力データ（p）を記憶する流動データ記憶部と、
圧力項を除いた流体力学における運動量保存の式を時間積分して得られる以下の式（１）を空間的に離散化し、前記流動場データ記憶部から読み出した流速データ（v）を代入して暫定的な流速（v^*）を計算する暫定流速計算部と、

前記運動量保存の式と前記質量保存の式とを連成し、かつ右辺第１項の前記暫定流速計算部により計算された暫定的な流速（v^*）に係る係数部にダルシー則抵抗を加えて得られる以下の式（２）を空間的に離散化し、前記流動場データ記憶部から読み出した圧力データ（p）および前記暫定流速計算部により計算された暫定的な流速（v^*）を代入して、圧力（p^*）が収束するまで繰り返し計算をする圧力計算部と、

以下の式（３）に、前記暫定流速計算部により計算された前記暫定的な流速（v^*）と前記圧力計算部により計算された圧力の値（p^*）とを代入して時間刻み幅（Δt）後の流速（v^k+1）に補正する流速補正部と

してコンピュータを機能させる。

また、本発明に係る電場解析プログラムは、電気化学反応により電力を発生させる陰極側多孔質層、陰極側触媒層、電解質膜、陽極側触媒層および陽極側多孔質層からなる電場の電子ポテンシャルおよびイオンポテンシャルを計算する電場解析プログラムであって、
電子ポテンシャル（φ_s）、イオンポテンシャル（φ_f）、前記電子ポテンシャル（φ_s）と前記イオンポテンシャル（φ_f）との差であるポテンシャル差（η）、陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を記憶する電場データ記憶部と、
以下の式（４）および式（５）に、前記電場データ記憶部から読み出した前記ポテンシャル差（η）を代入して暫定的な陰極側の輸送電流（Ｊ⁻）および暫定的な陽極側の輸送電流（Ｊ^＋）を計算する暫定輸送電流計算部と、

以下の式（６）を空間的に離散化し、外部電位を境界条件として、前記暫定的な陰極側の輸送電流（Ｊ⁻）および前記暫定的な陽極側の輸送電流（Ｊ^＋）を代入して暫定的な電子ポテンシャル（φ_s）を計算する暫定電子ポテンシャル計算部と、

以下の式（７）を空間的に離散化し、陰極側の境界条件を自由端とするとともに、陽極側の境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端として、前記暫定的な陰極側の輸送電流（Ｊ⁻）および前記暫定的な陽極側の輸送電流（Ｊ^＋）を代入して暫定的なイオンポテンシャル（φ_f）を計算する暫定イオンポテンシャル計算部と、

前記式（７）の逆関数に前記暫定的なイオンポテンシャル（φ_f）を代入してポテンシャル差（η）を計算するポテンシャル差計算部と、
前記式（４）および前記式（５）に前記ポテンシャル差計算部により得られたポテンシャル差（η）を代入して陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する輸送電流計算部と してコンピュータを機能させる。

さらに、本発明に係る電場解析プログラムは、電子ポテンシャル（φ_s）、イオンポテンシャル（φ_f）、前記電子ポテンシャル（φ_s）と前記イオンポテンシャル（φ_f）との差であるポテンシャル差（η）、陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を記憶する電場データ記憶部と、
前記式（４）および前記式（５）に、前記電場データ記憶部から読み出した前記ポテンシャル差（η）を代入して陰極側の固体層と触媒層との陰極側境界面における暫定的な輸送電流（Ｊ⁻）、および陽極側の固体層と触媒層との陽極側境界面における暫定的な輸送電流（Ｊ^＋）を計算する境界面暫定輸送電流計算部と、
以下の式（８）を空間的に離散化し、外部電位を境界条件として、前記陰極側境界面における輸送電流（Ｊ⁻）および陽極側境界面における輸送電流（Ｊ^＋）を代入して暫定的な電子ポテンシャル（φ_s）を計算する暫定電子ポテンシャル計算部と、

以下の式（９）を空間的に離散化し、陰極側境界面における境界条件を自由端とするとともに、陽極側の境界面における境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端として、前記陰極側境界面における輸送電流（Ｊ⁻）および陽極側境界面における輸送電流（Ｊ^＋）を代入して各境界面における暫定的なイオンポテンシャル（φ_f）を計算する暫定イオンポテンシャル計算部と、

前記式（９）の逆関数に陰極側境界面における前記暫定的なイオンポテンシャル（φ_f）を代入してポテンシャル差（η）を計算するポテンシャル差計算部と、
前記式（４）および前記式（５）に前記ポテンシャル差計算部により得られたポテンシャル差（η）を代入して陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する輸送電流計算部としてコンピュータを機能させる。

本発明によれば、流動場や電場におけるシミュレーション解析の解析速度を向上させたり、解析負荷を軽減させたりすることができる。

本発明に係る第１実施形態の流動場解析プログラムおよび電場解析プログラムを機能させるコンピュータ解析システムの構成を示すブロック図である。 燃料電池の概略構成を示す概略構成図である。 本第１実施形態における流動場解析機能の処理を示すフローチャートである。 本第１実施形態における電場解析機能の処理を示すフローチャートである。 本発明に係る第２実施形態の流動場解析プログラムおよび電場解析プログラムを機能させるコンピュータ解析システムの構成を示すブロック図である。 本第２実施形態における電場解析機能の処理を示すフローチャートである。 従来の流動場解析プログラムによる圧力計算の解析結果を示すグラフである。 本発明に係る第１実施形態の流動場解析プログラムによる圧力計算の解析結果を示すグラフである。 繰り返し計算の回数に対する従来の流動場解析プログラムと本発明に係る第１実施形態の流動場解析プログラムとによる圧力計算の解析結果の二乗平均の値を比較して示すグラフである。 従来の電場解析プログラムによる陽極側多孔質層内および陽極側触媒層内における電子ポテンシャルの解析結果を示すグラフである。 従来の電場解析プログラムによる陰極側多孔質層内および陰極側触媒層内における電子ポテンシャルの解析結果を示すグラフである。 従来の電場解析プログラムによる電解質膜内におけるイオンポテンシャルの解析結果を示すグラフである。 本発明に係る第１実施形態の電場解析プログラムによる陽極側多孔質層内および陽極側触媒層内における電子ポテンシャルの解析結果を示すグラフである。 本発明に係る第１実施形態の電場解析プログラムによる陰極側多孔質層内および陰極側触媒層内における電子ポテンシャルの解析結果を示すグラフである。 本発明に係る第１実施形態の電場解析プログラムによる電解質膜内におけるイオンポテンシャルの解析結果を示すグラフである。 繰り返し計算の回数に対する、従来の電場解析プログラムと本発明に係る第１実施形態の電場解析プログラムとによる陽極側多孔質層内および陽極側触媒層内における電子ポテンシャルの解析結果の二乗平均の値を比較して示すグラフである。 繰り返し計算の回数に対する、従来の電場解析プログラムと本発明に係る第１実施形態の電場解析プログラムとによる陰極側多孔質層内および陰極側触媒層内における電子ポテンシャルの解析結果の二乗平均の値を比較して示すグラフである。 繰り返し計算の回数に対する、従来の電場解析プログラムと本発明に係る第１実施形態の電場解析プログラムとによる電解質膜内におけるイオンポテンシャルの解析結果の二乗平均の値を比較して示すグラフである。 本発明に係る第２実施形態の電場解析プログラムによる陽極側多孔質層内における電子ポテンシャルの解析結果を示すグラフである。 本発明に係る第２実施形態の電場解析プログラムによる陰極側多孔質層内における電子ポテンシャルの解析結果を示すグラフである。 本発明に係る第２実施形態の電場解析プログラムによる電解質膜内におけるイオンポテンシャルの解析結果を示すグラフである。 繰り返し計算の回数に対する、本発明に係る第２実施形態の電場解析プログラムによる陽極側多孔質層内における電子ポテンシャルの解析結果の二乗平均の値を示すグラフである。 繰り返し計算の回数に対する、本発明に係る第２実施形態の電場解析プログラムによる陰極側多孔質層内における電子ポテンシャルの解析結果の二乗平均の値を示すグラフである。 繰り返し計算の回数に対する、本発明に係る第２実施形態の電場解析プログラムによる電解質膜内におけるイオンポテンシャルの解析結果の二乗平均の値を示すグラフである。

以下、本発明に係る流動場解析プログラムおよび電場解析プログラムの実施形態について図面を用いて説明する。

図１は、第１実施形態の流動場解析プログラム１ａおよび電場解析プログラム１ｂを機能させるコンピュータ解析システム１Ａの構成を示すブロック図である。図１に示すように、本第１実施形態のコンピュータ解析システム１Ａは、表示手段２、入力手段３、記憶手段４および演算処理手段５から構成される。

表示手段２は、液晶ディスプレイやCRTディスプレイ等から構成されており、流動場や電場に関する各種のデータや燃料電池内の流動や起電力のシミュレーション結果を表示するものである。また、入力手段３としては、マウスやキーボード等から構成されている。

記憶手段４は、各種のデータやプログラムを記憶するものであり、ハードディスク、フラッシュメモリおよびRAM（Random Access Memory）等から構成されている。本第１実施形態において、記憶手段４は、プログラム記憶部４１、流動場データ記憶部４２、電場データ記憶部４３、エネルギーデータ記憶部４４および化学種保存データ記憶部４５を有している。前記プログラム記憶部４１には、流動場解析プログラム１ａ、電場解析プログラム１ｂおよびオプション計算プログラム１ｃがインストールされている。また、前記流動場データ記憶部４２には、初期値としての流速データ（v）および圧力データ（p）が記憶されており、前記電場データ記憶部４３には、初期値としての電子ポテンシャル（φ_s）とイオンポテンシャル（φ_f）との差であるポテンシャル差（η）および外部電位の各種データが記憶されている。

演算処理手段５は、CPU（Central Processing Unit）等から構成されており、記憶手段４にインストールされた流動場解析プログラム１ａ、電場解析プログラム１ｂおよびオプション計算プログラム１ｃを実行させることにより、流動場解析機能６、電場解析機能７Ａおよびオプション機能８とを実現するようになっている。以下、各機能を実現するための各構成部についてより詳細に説明する。

流動場解析機能６は、以下の式（１０）に示す流体力学における質量保存の式と、式（１１）に示す流体力学における運動量保存の式とを連成させて、所定の流動場の流速や圧力を計算する機能であり、暫定流速計算部６１と、圧力計算部６２と、流速補正部６３とを有している。

なお、上記式（１１）は、右辺第１項が圧力項、右辺第２項が対流項、右辺第３項が外力項および右辺第４項がダルシー則抵抗を表す項である。

暫定流速計算部６１は、式（１１）の圧力項を除いて時間積分することにより得られる以下の式（１）を空間的に離散化した式を用い、前記流動場データ記憶部４２から読み出した初期値としての流速データ（v）を代入して暫定的な流速（v^*）を計算するものである。

圧力計算部６２は、前記式（１１）と前記式（１０）とを連成し、かつ右辺第１項の前記暫定流速計算部６１により計算された暫定的な流速（v^*）に係る係数部にダルシー則抵抗を加えて得られる以下の式（２）を空間的に離散化した式を用い、前記流動場データ記憶部４２から読み出した初期値としての圧力データ（p）および前記暫定流速計算部６１により計算された暫定的な流速（v^*）を代入して、圧力（p^*）が収束するまで繰り返し計算をするものである。

なお、従来の解析方法では、圧力や流速の時間変化が主に対流項に支配されていると仮定されている。しかし、多孔質材内の流路が狭い場所を流れる場合には、対流項が小さくダルシー則抵抗の項が支配的となり、この仮定は成り立たなくなる。そこで、本第１実施形態における圧力計算部６２では、圧力とダルシー則抵抗が力学的に釣り合う状況を仮定して、右辺第１項の流速（v^*）に係る係数部にダルシー則抵抗を直接的に加えることにより、圧力（p^*）の計算を行っている。なお、流動場に多孔質材が含まれない場合において、本第１実施形態における前記式（２）は、右辺第１項のダルシー則抵抗の分母の浸透率κが無限大、分子の多孔率εがゼロとなり、その値をゼロとして扱えるため、従来の解析方法の圧力計算で用いられる式と同じ式として扱える。

流速補正部６３は、以下の式（３）に、前記暫定流速計算部６１により計算された暫定的な流速（v^*）と前記圧力計算部６２により計算された圧力の値（p^*）とを代入して時間刻み幅（Δt）後の流速（v^k+1）に補正するものである。

なお、本第１実施形態における流動解析機能６では、前記流動場データ記憶部４２に記憶されている流速データ（v）および圧力データ（p）を、前記圧力計算部６２および前記流速補正部６３により計算された流速（v^k+1）および圧力（p^*）のデータに置き換え、それを初期値として、時間刻み幅（Δt）毎の解析が行われる。

つぎに、電場解析機能７Ａおよびその各構成部について詳細に説明する。

電場解析機能７Ａは、図２に示すような燃料電池における陰極側の多孔質層、この陰極側の多孔質層と電解質膜との間の陰極側触媒層、電解質膜、陽極側の多孔質層と電解質膜との間の陽極側触媒層および陽極側多孔質層からなる電場において、電子ポテンシャル（φ_s）、イオンポテンシャル（φ_f）、ポテンシャル差（η）、陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する機能である。

本第１実施形態の電場解析機能７Ａは、暫定輸送電流計算部７１と、暫定電子ポテンシャル計算部７２と、暫定イオンポテンシャル計算部７３と、ポテンシャル差計算部７４と、輸送電流計算部７５とを有している。

暫定輸送電流計算部７１は、陰極側の輸送電流（Ｊ⁻）を示す以下の式（４）と、陽極側の輸送電流（Ｊ^＋）を示す以下の式（５）とに、前記電場データ記憶部４３から読み出した初期値としての前記ポテンシャル差（η）を代入して暫定的な陰極側の輸送電流（Ｊ⁻）および暫定的な陽極側の輸送電流（Ｊ^＋）を計算するものである。

なお、図２に示す陰極側流路および陰極側多孔質層を流れる気体が水素であり、陽極側流路および陽極側多孔質層を流れる気体が酸素である場合、上記式（４）および上記式（５）の各基準電流（i₀ ^±）は、以下の式（１２）および式（１３）により算出される定数値である。

また、前記式（５）の自由エネルギーポテンシャル（U₀）は以下の式（１４）により算出される定数値である。

暫定電子ポテンシャル計算部７２は、以下の式（６）を空間的に離散化した式に、境界条件として外部電位を与え、かつ前記暫定的な陰極側の輸送電流（Ｊ⁻）および前記暫定的な陽極側の輸送電流（Ｊ^＋）を代入して暫定的な電子ポテンシャル（φ_s）を計算するものである。この式（６）は、陽極側触媒層および陰極側触媒層において、電子ポテンシャル（φ_s）の空間的な勾配と各輸送電流（Ｊ^±）とが等しくなることを示すものであり、陽極側触媒層では前記式（５）で求められる輸送電流（Ｊ^＋）が発生し、陰極側触媒層では前記式（４）で求められる輸送電流（Ｊ⁻）が消滅し、その他の領域では電子の発生や消滅が無いため電子ポテンシャル（φ_s）の勾配がゼロであることを示している。なお、外部電位は、図２に示す燃料電池により起動される外部機器に必要とされる電位であり、任意に設定されるものである。

暫定イオンポテンシャル計算部７３は、以下の式（７）を空間的に離散化した式に、陰極側の境界条件を自由端として与えるとともに、陽極側の境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端として与え、かつ前記暫定的な陰極側の輸送電流（Ｊ⁻）および前記暫定的な陽極側の輸送電流（Ｊ^＋）を代入して暫定的なイオンポテンシャル（φ_f）を計算するものである。この式（７）は、陽極側触媒層および陰極側触媒層において、イオンポテンシャル（φ_f）の空間的な勾配と各輸送電流（Ｊ^±）とが等しくなることを示すものであり、陽極側触媒層では式（５）で求められる輸送電流（Ｊ^＋）が消滅し、陰極側触媒層では前記式（４）求められる輸送電流（Ｊ⁻）が発生し、その他の解析領域ではイオンの発生や消滅が無いためイオンポテンシャル（φ_f）の勾配がゼロであることを示している。

なお、従来の解析方法では、このイオンポテンシャル（φ_f）の計算において、陰極側と陽極側の両側の境界条件を自由端として計算をしている。しかし、両側の境界条件を自由端とすると、得られるイオンポテンシャル（φ_f）は暫定電子ポテンシャル計算部で代入される初期値としてのポテンシャル差（η）に依存してしまう。そのため、得られるイオンポテンシャル（φ_f）が正しいか否かは、実験値との比較等を行わなければならない。そこで、本第１実施形態では、陽極側の境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端として与えている。

ポテンシャル差計算部７４は、前記式（４）と前記式（７）とを連成させた式について、さらにポテンシャル差（η）が解となる逆関数の式を求め、前記暫定的なイオンポテンシャル（φ_f）を代入してポテンシャル差（η）を計算する。なお、前記式（５）と前記式（７）とを連成させて、その逆関数の式からもポテンシャル差（η）を計算することも可能であるが、前記式（５）は指数関数で表現されており、逆関数を求める際に数値的な誤差が含まれてしまう。よって、計算精度の点から、本第１実施形態のように前記式（４）を用いることが望ましい。

輸送電流計算部７５は、前記式（４）および前記式（５）に前記ポテンシャル差計算部により得られたポテンシャル差（η）を代入して陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算するものである。

つぎに、オプション機能８について説明する。オプション機能８とは、前記流動場解析機能６や前記電場解析機能７Ａにより得られた、圧力、流速、電子ポテンシャル、イオンポテンシャル等に基づいて、各種物理量を計算するための機能である。本第１実施形態において、オプション機能８は、エネルギー計算部８１と、化学種保存式計算部８２とを有している。

エネルギー計算部８１は、多孔質層内のジュール熱や電気化学反応による仕事量のエネルギーバランスを計算するためのものである。本第１実施形態におけるエネルギー計算部８１は、以下の式（１５）に前記流動場解析機能６により得られた流速（v）および圧力（p）と、前記電場解析機能７Ａにより得られた各輸送電流（Ｊ^±）およびポテンシャル差（η）とを代入して所定の領域のエネルギーを計算することが可能である。

化学種保存式計算部８２は、各触媒層における化学反応の生成を計算するためのものであり、以下の式（１６）に前記流動場解析機能６により得られた流速（v）と、前記電場解析機能７Ａにより得られた各輸送電流（Ｊ^±）とを代入して化学種の保存を計算することが可能である。

なお、オプション機能８については、計算させる内容が上記エネルギー計算等に限定されるものではなく、その他の温度場の計算等をさせてもよい。

つぎに、本第１実施形態の流動場解析機能６、電場解析機能７Ａおよびオプション機能８の作用について、図３および図４のフローチャートを参照しつつ説明する。

まず、図３を参照しつつ流動場解析機能６の処理フローについて説明する。ステップＳ１１では、暫定流速計算部６１が、流動場データ記憶部４２から流速データ（v）を読み出し、この流速データ（v）を式（１）に代入して流速（v^*）を計算し、これを暫定的な流速（v^*）として流動場データ記憶部４２に記憶させる。

つぎに、ステップＳ１２では、圧力計算部６２が、流動場データ記憶部４２から圧力データ（p）を読み出し、かつ暫定流速計算部６１によりステップＳ１１で得られた暫定的な流速（v^*）を式（２）に代入して圧力（p^*）を計算する。

ステップＳ１３では、ステップＳ１２の処理における圧力（p^*）の計算が収束したか否かを判別する。本第１実施形態においては、流動場データ記憶部４２から読み出された圧力データ（p）とステップＳ１２で計算された圧力（p^*）との差の合計が所定の値より小さいか否かで判別する。当該圧力の差の合計が所定の値より小さい場合には、圧力（p^*）の計算が収束したと判別し、計算された圧力（p^*）を前記流動場データ記憶部４２の圧力データ（p）に置き換えてステップＳ１４へ進む（ステップＳ１３：ＹＥＳ）。一方、当該圧力の差の合計が所定の値より大きい場合には、圧力（p^*）の計算が収束していないと判別し、計算された圧力（p^*）を前記流動場データ記憶部４２の圧力データ（p）に置き換えてステップＳ１２へ戻り（ステップＳ１３：ＮＯ）、収束するまで圧力（p^*）を求める計算と収束判別処理が繰り返される。

なお、本第１実施形態における式（２）には、流速（v^*）に係る係数部にダルシー則抵抗の項が直接的に加えられているため、対流項が小さくダルシー則抵抗の項が支配的な場合においてもダルシー則抵抗の影響を直接的に圧力計算に反映させることとなり、ステップＳ１２〜ステップＳ１３の収束性が向上する。また、前述のように、式（２）に加えられたダルシー則抵抗の項は、対流項が支配的な場合にはゼロとして扱えるため、従来の方法の計算方法と等価になる。

ステップＳ１４では、暫定流速計算部６１により計算された暫定的な流速（v^*）と圧力計算部６２により計算された圧力の値（p^*）とを代入して時間刻み幅（Δt）後の流速（v^k+1）に補正する計算を行う。以上の各ステップＳ１１〜Ｓ１４により時間刻み幅（Δt）後の流速（v^k+1）および圧力の値（p^*）を計算することができる。

なお、所定時間経過後の流速（v）および圧力（p）を得るには、ステップＳ１４の後、前記流動場データ記憶部４２の流速データ（v）および圧力データ（p）を、前記流速補正部６３および前記圧力計算部６２により計算された流速（v^k+1）および圧力（p^*）のデータに置き換え、各ステップＳ１１〜Ｓ１４を所定時間まで順次計算することになる。

つぎに、電場解析機能７Ａの各構成部における作用について図４のフローチャートを参照しつつ説明する。まずステップＳ２１において、暫定輸送電流計算部７１が、電場データ記憶部４３からポテンシャル差（η）を読み出し、このポテンシャル差（η）を式（４）および式（５）に代入して、暫定的な陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する。

つぎにステップＳ２２では、暫定電子ポテンシャル計算部７２が、電場データ記憶部４３から外部電位を読み出し、この外部電位、前記ステップＳ２１で計算された暫定的な陰極側の輸送電流（Ｊ⁻）および暫定的な陽極側の輸送電流（Ｊ^＋）をそれぞれ式（６）に代入して暫定的な電子ポテンシャル（φ_s）を計算する。

また、ステップＳ２３では、暫定イオンポテンシャル計算部７３が、式（７）に暫定的な陰極側の輸送電流（Ｊ⁻）および暫定的な陽極側の輸送電流（Ｊ^＋）を代入して暫定的なイオンポテンシャル（φ_f）を計算する。このとき、陰極側の境界条件を自由端とするとともに、陽極側の境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端とする。

ステップＳ２４では、ポテンシャル差計算部７４が、前記ステップ２３で得られた暫定的なイオンポテンシャル（φ_f）を式（４）と式（７）とを連成させた式の逆関数の式に代入してポテンシャル差（η）を計算する。

ステップＳ２５では、輸送電流計算部７５が、ステップＳ２４で得られたポテンシャル差（η）を前記式（４）および前記式（５）に代入することにより陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する。

なお、本第１実施形態における電場解析機能７Ａでは、ステップＳ２６において、ステップ２４の処理におけるポテンシャル差（η）の計算が収束したか否かを判別する。本第１実施形態においては、電場データ記憶部から読み出されたポテンシャル差（η）とステップ２４で計算されたポテンシャル差（η）との差の合計が所定の値より小さいか否かで判別する。当該ポテンシャル差（η）の差の合計が所定の値より小さい場合には、ポテンシャル差（η）の計算が収束したと判別し、計算されたポテンシャル差（η）を前記電場データ記憶部４３のポテンシャル差（η）に置き換えて計算を終了する（ステップＳ２６：ＹＥＳ）。一方、当該ポテンシャル差（η）の差の合計が所定の値より大きい場合には、ポテンシャル差（η）の計算が収束していないと判別し、計算されたポテンシャル差（η）を前記電場データ記憶部４３のポテンシャル差（η）に置き換えてステップＳ２１へ戻り（ステップＳ２６：ＮＯ）、収束するまで各種データを求める計算と収束判別処理が繰り返される。

つぎに、オプション機能８の各構成部における作用について説明する。エネルギー計算部８１は、上記ステップＳ１１〜ステップＳ１４で求めた流動場における流速および圧力と、上記ステップＳ２１〜ステップＳ２６で求めた電場におけるポテンシャル差（η）、陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）とを式（１５）に代入し、流動場と電場とが等しく存在する場合のエネルギーを計算する。

また、化学種保存式計算部８２は、上記ステップＳ１１〜ステップＳ１４で求めた流動場における流速と、上記ステップＳ２１〜ステップＳ２６で求めた電場における陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）とを式（１６）に代入し、流動場と電場とが等しく存在する場合の化学種の保存に関する値を計算する。

なお、本第１実施形態において、流動場解析機能６、電場解析機能７Ａと順次に行った後にエネルギー計算８１等を行ったが、計算順はこれに限定されるものではなく、例えば、流動場解析機能６と電場解析機能７Ａの解析順は逆でもよく、流動場解析機能６と電場解析機能７Ａとを同時並列に計算するようにしてもよい。

つぎに、本発明に係る電場解析プログラムの第２実施形態について説明する。なお、本第２実施形態のうち、前述した第１実施形態の構成と同一若しくは相当する構成については同一の符号を付して再度の説明を省略する。

図５は、本第２実施形態における、流動場解析プログラム１ａおよび電場解析プログラム１ｄとして機能するコンピュータ解析システム１Ｂの構成を示すブロック図である。本第２実施形態のコンピュータ解析システム１Ｂは、表示手段２、入力手段３、記憶手段４および演算処理手段５から構成されている。本第２実施形態の特徴は、前述した第１実施形態における演算処理手段５の電場解析機能７Ｂにおいて、暫定輸送電流計算部７１に代えて、触媒層と多孔質層との境界における暫定的な輸送電流を算出するための境界面暫定輸送電流計算部７６を備えた点にある。

すなわち、本第２実施形態における電場解析機能７Ｂは、陰極側触媒層および陽極側触媒層が十分に薄く、かつ電気ポテンシャル（φ_s）が陰極側触媒層と陰極側多孔質層の境界および陽極側触媒層と陽極側多孔質層の境界で不連続に変化すると仮定することにより、電場の格子間隔を各触媒層と等しいかそれ以上の幅に設定し、各触媒層内の計算をせずに電子ポテンシャル（φ_s）、イオンポテンシャル（φ_f）、ポテンシャル差（η）、陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算するためのものである。このため、本第２実施形態における電場解析機能７Ｂは、境界面暫定輸送電流計算部７６と、暫定電子ポテンシャル計算部７２と、暫定イオンポテンシャル計算部７３と、ポテンシャル差計算部７４と、輸送電流計算部７５とを有している。

境界面暫定輸送電流計算部７６は、前記式（４）および前記式（５）に、前記電場データ記憶部４３から読み出した初期値としての前記ポテンシャル差（η）を代入して陰極側の固体層と触媒層との陰極側境界面における暫定的な輸送電流（Ｊ⁻）、および陽極側の固体層と触媒層との陽極側境界面における暫定的な輸送電流（Ｊ^＋）を計算するものである。

本第２実施形態における暫定電子ポテンシャル計算部７２は、以下の式（８）を空間的に離散化した式に、境界条件として外部電力を与え、前記境界面暫定輸送電流計算部７６により得られた陰極側境界面における暫定的な輸送電流（Ｊ⁻）および陽極側境界面における暫定的な輸送電流（Ｊ^＋）を代入して暫定的な電子ポテンシャル（φ_s）を計算するものである。

本第２実施形態における暫定イオンポテンシャル計算部７３は、以下の式（９）を空間的に離散化した式に、陰極側境界面における境界条件を自由端とするとともに、陽極側の境界面における境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端として、前記陰極側境界面における暫定的な輸送電流（Ｊ⁻）および陽極側境界面における暫定的な輸送電流（Ｊ^＋）を代入して各境界面における暫定的なイオンポテンシャル（φ_f）を計算するものである。

すなわち、本第２実施形態の電場解析機能７Ｂは、電子ポテンシャル（φ_s）およびイオンポテンシャル（φ_f）の発生および消滅が陰極側境界面および陽極側境界面で不連続であとして、触媒層内の解を求めないことにより、計算負荷を軽減したものである。

本第２実施形態におけるポテンシャル差計算部７４は、前記式（４）と前記式（９）とを連成させ、ポテンシャル差（η）が解となる逆関数に、陰極側境界面における前記暫定的なイオンポテンシャル（φ_f）を代入してポテンシャル差（η）を計算するものである。

本第２実施形態における輸送電流計算部７５は、前記式（４）および前記式（５）に本第２実施形態におけるポテンシャル差計算部７４により得られたポテンシャル差（η）を代入して陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算するものである。

つぎに、本第２実施形態における電場解析機能７Ｂの各構成部における作用について図６のフローチャートを参照しつつ説明する。

ステップＳ３１では、境界面暫定輸送電流計算部７６が、電場データ記憶部４３からポテンシャル差（η）を読み出し、このポテンシャル差（η）を式（４）および式（５）に代入して、陰極側境界面における暫定的な輸送電流（Ｊ⁻）および暫定的な陽極側境界面における輸送電流（Ｊ^＋）を計算する。

つぎにステップＳ３２では、暫定電子ポテンシャル計算部７２が、電場データ記憶部４３から外部電位を読み出し、この外部電位と、前記ステップＳ３１で計算された陰極側境界面における暫定的な輸送電流（Ｊ⁻）および陽極側境界面における暫定的な輸送電流（Ｊ^＋）を式（８）に代入して暫定的な電子ポテンシャル（φ_s）を計算する。

ステップＳ３３では、暫定イオンポテンシャル計算部７３が、式（９）に前記陰極側境界面における暫定的な輸送電流（Ｊ⁻）および前記陽極側境界面における暫定的な輸送電流（Ｊ^＋）を代入して暫定的なイオンポテンシャル（φ_f）を計算する。ただし、陰極側境界面における境界条件を自由端とするとともに、陽極側の境界面における境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端とする。

ステップＳ３４では、ポテンシャル差計算部７４が、前記ステップ３３で得られた暫定的なイオンポテンシャル（φ_f）を式（４）と式（９）とを連成させた式の逆関数の式に代入してポテンシャル差（η）を計算する。

ステップＳ３５では、輸送電流計算部７５が、ステップＳ３４で得られたポテンシャル差（η）を、前記式（４）および前記式（５）に代入することにより陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する。

なお、本第２実施形態における電場解析機能７Ｂでは、ステップＳ３６において、ステップ３４の処理におけるポテンシャル差（η）の計算が収束したか否かを判別する。本第２実施形態においては、電場データ記憶部から読み出されたポテンシャル差（η）とステップ３４で計算されたポテンシャル差（η）との差の合計が所定の値より小さいか否かで判別する。当該ポテンシャル差（η）の差の合計が所定の値より小さい場合には、ポテンシャル差（η）の計算が収束したと判別し、計算されたポテンシャル差（η）を前記電場データ記憶部４３のポテンシャル差（η）に置き換えて計算を終了する（ステップＳ３６：ＹＥＳ）。一方、当該ポテンシャル差（η）の差の合計が所定の値より大きい場合には、ポテンシャル差（η）の計算が収束していないと判別し、計算されたポテンシャル差（η）を前記電場データ記憶部４３のポテンシャル差（η）に置き換えてステップＳ３１へ戻り（ステップＳ３６：ＮＯ）、収束するまで各種データを求める計算と収束判別処理が繰り返される。

以上のような本第１実施形態および第２実施形態によれば、以下の効果が得られる。
１．多孔質材を含む流動場における計算の収束性を向上することができる。
２．電場の解析において、初期値に対する解の依存性を軽減することができる。
３．電場の解析において、格子解像度をある程度大きくとることにより、計算負荷を軽減することができる。
このような作用効果を奏することにより、燃料電池などの電気化学反応装置における電気化学ポテンシャルと流動が連成する問題において収束性がよく、計算速度の速いシミュレーションが実現できる。

つぎに、本発明に係る流動場解析プログラムの実施例を図面を用いて説明する。

本実施例１では、従来の方法の流動解析プログラムと、本発明に係る流動解析プログラムとにより多孔質材を含む流動場の計算を行い、圧力計算の収束性を比較検討した。ただし、流体の密度（ρ）を1.2kg/m³、粘性係数（μ）を1.0*10^-5Pa・s、多孔質材の多孔率（ε）を0.8、多孔質材の浸透率（κ）を1.0*10^-11m²、出口の圧力データ（p）を0N/m²、入口の流速データ（v）を0.005m/s、時間刻み幅（Δt）を1.0*10^-6s、格子間隔（Δx）を1.0*10^-5m、格子数（N）を1000とし、格子数（N）が325〜750の間を多孔質材として計算を行った。

図７および図８は、従来の方法の流動解析プログラムと、本発明に係る流動解析プログラムとにより、それぞれ100000回の繰り返し計算を行ったときの圧力分布を示すグラフであって、繰り返し数が1000回目、5000回目、10000回目、15000回目、20000回目、25000回目、30000回目および100000回目の圧力分布を示している。なお、図７および図８の横軸は格子数（N）を示しており、縦軸は圧力（p^*）を示している。本実施例１では、従来の方法の流動解析プログラムおよび本発明に係る流動解析プログラムともに、100000回繰り返し計算した結果を十分に収束したものとした。

図７および図８に示すように、従来の方法の流動解析プログラムの解析結果方および本発明に係る流動解析プログラムの解析結果のいずれも、繰り返し計算の回数が増加するにつれて収束値に漸近していく。しかし、本発明の流動解析結果の方がより小さい回数で収束値に近づいており、たとえば繰り返し計算の回数30000回目を比較すると、本発明の流動解析結果は収束値と略同一の値まで漸近しているのに対し、従来法による流動解析結果は収束値との差が存在し、収束していない。

ここで、従来法の流動解析プログラムと、本発明に係る流動解析プログラムとの収束性の差を、収束値と各結果の値との差の二乗平均により示す。

図９は、従来法の流動解析結果と、本発明の流動解析結果とにおける前記二乗平均の値を示している。横軸は繰り返し計算の回数を対数で示しており、縦軸は前記二乗平均を対数で示している。図９に示すように、従来法の流動解析結果および本発明の流動解析結果のいずれも、繰り返し計算の回数が増加するにつれて二乗平均の値が次第に小さくなっている。しかし、本発明の流動解析結果の方が従来法の流動解析結果に比べて前記二乗平均の値がより早く小さくなっている。例えば、前記二乗平均の値が0.0001以下になるための繰り返し計算の回数は、本発明の流動解析結果では、10000回以下であるのに対し、従来法の流動解析結果では、約30000回以上でなければならないことが確認された。

以上に示したように、本実施例１における図７〜図９に示した結果によると、多孔質材を含む流動場の計算において、本発明に係る流動解析プログラムは、従来の方法の流動解析プログラムに比べ、収束性が高いことが示された。

つぎに、本発明に係る電場解析プログラムの実施例２を図面を用いて説明する。本実施例２では、従来の方法の電場解析プログラムと、本発明に係る第１実施形態の電場解析プログラムと、本発明に係る第２実施形態の電場解析プログラムとにより陽極側多孔質層、陽極側触媒層、電解質膜、陰極側触媒層および陰極側多孔質層からなる電場の計算を行い、各多孔質層内および各触媒層内の電子ポテンシャル（φ_s）と、電解質膜内のイオンポテンシャル（φ_f）の計算の収束性を比較検討した。

ただし、各プログラムにおいて、電子コンダクティビティ（σ_s）を300Ω^-1m^-1、イオンコンダクティビティ（σ_f）を10Ω^-1m^-1、外部電位を0.65V、アノード移動係数（α_a）を10⁸A/m³、カソード移動係数（α_c）を10⁴A/m³、各多孔質層の厚さを50*10^-6m、触媒層の厚さ（δ）を5*10^-6m、電解質膜の厚さを50*10^-6m、格子間隔（Δx）を0.5*10^-6m、温度を353.0Kとして計算を行った。

図１０は、従来の方法による電場解析プログラムの解析結果の陽極側多孔質層および陽極側触媒層の電子ポテンシャル（φ_s）を示すものである。また、図１１は、従来法による電場解析結果の陰極側多孔質層および陰極側触媒層の電子ポテンシャル（φ_s）を示すものであり、図１２は従来法による電場解析結果の電解質膜のイオンポテンシャル（φ_f）を示すものである。なお、図１０〜図１２の横軸は格子数（N）を示しており、縦軸は、図１０および図１１では電子ポテンシャル（φ_s）、図１２ではイオンポテンシャル（φ_f）を示している。

図１０〜図１２に示すように、従来法による電場解析結果では、電子ポテンシャル（φ_s）およびイオンポテンシャル（φ_f）は、繰り返し計算の回数が10000000回目であっても十分に収束することはなかった。特に、電解質膜のイオンポテンシャル（φ_f）に関しては、収束せずに解が発散してしまうという結果になった。

以上のような従来法の電場解析結果に対し、図１３は、本発明に係る第１実施形態の電場解析プログラムによる解析結果において陽極側多孔質層および陽極側触媒層の電子ポテンシャル（φ_s）を示すものである。また、図１４は、第１実施形態による電場解析結果において陰極側多孔質層および陰極側触媒層の電子ポテンシャル（φ_s）を示すものであり、図１５は第１実施形態による電場解析結果において電解質膜のイオンポテンシャル（φ_f）を示すものである。

図１３〜図１５に示すように、第１実施形態による電場解析結果は、陽極側多孔質層および陽極側触媒層の電子ポテンシャル（φ_s）、陰極側多孔質層および陰極側触媒層の電子ポテンシャル（φ_s）および電解質膜のイオンポテンシャル（φ_f）のいずれの場合においても解の収束が確認された。

ここで、前記実施例１と同様に、従来法の電場解析プログラムと本発明に係る第１実施形態の電場解析プログラムとの収束性の差を、収束値と各結果の値との差の二乗平均により示す。

図１６は、従来法による電場解析結果および第１実施形態による電場解析結果における、陽極側多孔質層および陽極側触媒層の電子ポテンシャル（φ_s）の前記二乗平均の値である。なお、横軸は繰り返し計算の回数を対数で示しており、縦軸は前記二乗平均を対数で示している。同様に、図１７は各結果における陰極側多孔質層および陰極側触媒層の電子ポテンシャル（φ_s）の前記二乗平均の値、図１８は各結果における電解質膜のイオンポテンシャル（φ_f）の前記二乗平均の値を示している。なお、前記収束値には、10000000回繰り返し計算したときの第１実施形態による電場解析結果を用いた。

これら図１６〜図１８のいずれの場合においても見られるように、第１実施形態による電場解析結果の前記二乗平均の値が、従来法による電場解析結果の前記二乗平均の値に比べて明らかに早く小さくなっている。すなわち、本発明に係る第１実施形態の電場解析プログラムは、従来法の電場解析プログラムと比べて極めて収束性が高いことが確認された。

つぎに、本発明に係る第２実施形態の電場解析プログラムにより、陽極側触媒層内および陰極側触媒層内の計算を行わずに陽極側多孔質層の電子ポテンシャル（φ_s）、陰極側多孔質層の電子ポテンシャル（φ_s）および電解質膜のイオンポテンシャル（φ_f）の解析を行った。

図１９は、本発明に係る第２実施形態の電場解析プログラムによる解析結果において陽極側多孔質層の電子ポテンシャル（φ_s）を示すものである。また、図２０は、第２実施形態による電場解析結果において陰極側多孔質層の電子ポテンシャル（φ_s）を示すものであり、図２１は第２実施形態による電場解析結果において電解質膜のイオンポテンシャル（φ_f）を示すものである。

また、図２２は、第２実施形態による電場解析結果において陽極側多孔質層および陽極側触媒層の電子ポテンシャル（φ_s）の前記二乗平均の値であり、図２３は第２実施形態による電場解析結果において陰極側多孔質層の電子ポテンシャル（φ_s）の前記二乗平均の値、図２４は第２実施形態による電場解析結果において電解質膜のイオンポテンシャル（φ_f）の前記二乗平均の値を示している。なお、前記収束値には、10000000回繰り返し計算したときの解析結果を用いた。

図１９〜図２１に示すように、第２実施形態による電場解析結果は、第１実施形態の電場解析結果と同様に、陽極側多孔質層の電子ポテンシャル（φ_s）、陰極側多孔質層の電子ポテンシャル（φ_s）および電解質膜のイオンポテンシャル（φ_f）のいずれの場合においても、解析結果の収束が確認された。

また、図２１と図１８に示す電解質膜のイオンポテンシャル（φ_f）における二乗平均の値を例に、本発明に係る第１実施形態の電場解析プログラムおよび本発明に係る第２実施形態の電場解析プログラムの収束性について比較を行った。

それぞれの、繰り返し計算の回数が100000回目の二乗平均の値を比較すると、第１実施形態による電場解析結果の場合は二乗平均値が10^-6 程度であるのに対し、第２実施形態による電場解析結果の場合は二乗平均値が10^-7 程度であることがわかった。

すなわち、本実施例２においては、第２実施形態の電場解析プログラムの方が、第１実施形態の電場解析プログラムより収束性が高いことが確認された。

なお、本発明に係る流動場解析プログラムは、前述した実施形態に限定されるものではなく、適宜変更することができる。

例えば、式（２）に加えられるダルシー則抵抗の項は、高次精度で離散化したものを使用してもよい。

１Ａ、１Ｂ コンピュータ解析システム
１ａ 流動解析プログラム
１ｂ、１ｄ 電場解析プログラム
１ｃ オプション計算プログラム
２ 表示手段
３ 入力手段
４ 記憶手段
５ 演算処理手段
６ 流動場解析機能
７Ａ、７Ｂ 電場解析機能
８ オプション機能
４１ プログラム記憶部
４２ 流動場データ記憶部
４３ 電場データ記憶部
４４ エネルギーデータ記憶部
４５ 化学種保存データ記憶部
６１ 暫定流速計算部
６２ 圧力計算部
６３ 流速補正部
７１ 暫定輸送電流計算部
７２ 暫定電子ポテンシャル計算部
７３ 暫定イオンポテンシャル計算部
７４ ポテンシャル差計算部
７５ 輸送電流計算部
８１ エネルギー計算部
８２ 化学種保存式計算部

Claims (3)

1. 流体力学における運動量保存の式と流体力学における質量保存の式とを連成させて所定の流動場の流速および圧力を計算する流動場解析プログラムであって、
   流速データ（v）および圧力データ（p）を記憶する流動データ記憶部と、
   圧力項を除いた流体力学における運動量保存の式を時間積分して得られる以下の式（１）を空間的に離散化し、前記流動場データ記憶部から読み出した流速データ（v）を代入して暫定的な流速（v^*）を計算する暫定流速計算部と、
   

   

   前記運動量保存の式と前記質量保存の式とを連成し、かつ右辺第１項の前記暫定流速計算部により計算された暫定的な流速（v^*）に係る係数部にダルシー則抵抗を加えて得られる以下の式（２）を空間的に離散化し、前記流動場データ記憶部から読み出した圧力データ（p）および前記暫定流速計算部により計算された暫定的な流速（v^*）を代入して、圧力（p^*）が収束するまで繰り返し計算をする圧力計算部と、
   

   

   以下の式（３）に、前記暫定流速計算部により計算された前記暫定的な流速（v^*）と前記圧力計算部により計算された圧力の値（p^*）とを代入して時間刻み幅（Δt）後の流速（v^k+1）に補正する流速補正部と
   

   

   してコンピュータを機能させる流動場解析プログラム。
2. 電気化学反応により電力を発生させる陰極側多孔質層、陰極側触媒層、電解質膜、陽極側触媒層および陽極側多孔質層からなる電場の電子ポテンシャルおよびイオンポテンシャルを計算する電場解析プログラムであって、
   電子ポテンシャル（φ_s）、イオンポテンシャル（φ_f）、前記電子ポテンシャル（φ_s）と前記イオンポテンシャル（φ_f）との差であるポテンシャル差（η）、陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を記憶する電場データ記憶部と、
   以下の式（４）および式（５）に、前記電場データ記憶部から読み出した前記ポテンシャル差（η）を代入して暫定的な陰極側の輸送電流（Ｊ⁻）および暫定的な陽極側の輸送電流（Ｊ^＋）を計算する暫定輸送電流計算部と、
   

   

   

   以下の式（６）を空間的に離散化し、外部電位を境界条件として、前記暫定的な陰極側の輸送電流（Ｊ⁻）および前記暫定的な陽極側の輸送電流（Ｊ^＋）を代入して暫定的な電子ポテンシャル（φ_s）を計算する暫定電子ポテンシャル計算部と、
   

   

   以下の式（７）を空間的に離散化し、陰極側の境界条件を自由端とするとともに、陽極側の境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端として、前記暫定的な陰極側の輸送電流（Ｊ⁻）および前記暫定的な陽極側の輸送電流（Ｊ^＋）を代入して暫定的なイオンポテンシャル（φ_f）を計算する暫定イオンポテンシャル計算部と、
   

   

   前記式（７）の逆関数に前記暫定的なイオンポテンシャル（φ_f）を代入してポテンシャル差（η）を計算するポテンシャル差計算部と、
   前記式（４）および前記式（５）に前記ポテンシャル差計算部により得られたポテンシャル差（η）を代入して陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する輸送電流計算部と
   してコンピュータを機能させる電場解析プログラム。
3. 電気化学反応により電力を発生させる陰極側多孔質層、陰極側触媒層、電解質膜、陽極側触媒層および陽極側多孔質層からなる電場の電子ポテンシャルおよびイオンポテンシャルを計算する電場解析プログラムであって、
   電子ポテンシャル（φ_s）、イオンポテンシャル（φ_f）、前記電子ポテンシャル（φ_s）と前記イオンポテンシャル（φ_f）との差であるポテンシャル差（η）、陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を記憶する電場データ記憶部と、
   以下の式（４）および式（５）に、前記電場データ記憶部から読み出した前記ポテンシャル差（η）を代入して陰極側の固体層と触媒層との陰極側境界面における暫定的な輸送電流（Ｊ⁻）、および陽極側の固体層と触媒層との陽極側境界面における暫定的な輸送電流（Ｊ^＋）を計算する境界面暫定輸送電流計算部と、
   

   

   

   以下の式（８）を空間的に離散化し、外部電位を境界条件として、前記陰極側境界面における輸送電流（Ｊ⁻）および陽極側境界面における輸送電流（Ｊ^＋）を代入して暫定的な電子ポテンシャル（φ_s）を計算する暫定電子ポテンシャル計算部と、
   

   

   以下の式（９）を空間的に離散化し、陰極側境界面における境界条件を自由端とするとともに、陽極側の境界面における境界条件を前記暫定的な電子ポテンシャル（φ_s）から前記ポテンシャル差（η）を引いた固定端として、前記陰極側境界面における輸送電流（Ｊ⁻）および陽極側境界面における輸送電流（Ｊ^＋）を代入して各境界面における暫定的なイオンポテンシャル（φ_f）を計算する暫定イオンポテンシャル計算部と、
   

   

   前記式（９）の逆関数に陰極側境界面における前記暫定的なイオンポテンシャル（φ_f）を代入してポテンシャル差（η）を計算するポテンシャル差計算部と、
   前記式（４）および前記式（５）に前記ポテンシャル差計算部により得られたポテンシャル差（η）を代入して陰極側の輸送電流（Ｊ⁻）および陽極側の輸送電流（Ｊ^＋）を計算する輸送電流計算部と
   してコンピュータを機能させる電場解析プログラム。

Images