輸送方程式解析プログラムおよび輸送方程式解析方法
【課題】 解析精度の向上、解析速度の向上、解析負荷の軽減ならびに収束性の安定を図ることのできる輸送方程式解析プログラムおよび輸送方程式解析方法を提供する。
【解決手段】 空間的に離散化した各格子点における物理量を取得する物理量取得手段６と、第一領域側から境界へ流入出する物理量の流束と第二領域側から境界へ流入出する物理量の流束との連続性を示す流束境界条件式、および第一領域の輸送方程式に基づいて所定時間後の第一領域内における物理量の数値解を求める第一領域内解析手段７と、境界における第一領域側の物理量と第二領域側の物理量との連続性を示す物理量境界条件式、および第二領域の輸送方程式に基づいて所定時間後の第二領域内における物理量の数値解を求める第二領域内解析手段８とを有する。
【解決手段】 空間的に離散化した各格子点における物理量を取得する物理量取得手段６と、第一領域側から境界へ流入出する物理量の流束と第二領域側から境界へ流入出する物理量の流束との連続性を示す流束境界条件式、および第一領域の輸送方程式に基づいて所定時間後の第一領域内における物理量の数値解を求める第一領域内解析手段７と、境界における第一領域側の物理量と第二領域側の物理量との連続性を示す物理量境界条件式、および第二領域の輸送方程式に基づいて所定時間後の第二領域内における物理量の数値解を求める第二領域内解析手段８とを有する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、それぞれ性質の異なる複数の領域を移動する所定の物理量の輸送方程式を解析する輸送方程式解析プログラムおよび輸送方程式解析方法に関するものである。
【背景技術】
【０００２】
近年、燃料電池は、電力発生に伴うＣＯ２排出のないクリーンな次世代エネルギーとして注目されている。この燃料電池は、多孔質層や電解膜等から構成されており、水素等の燃料剤と酸素等の酸化剤を化学反応させることにより電力を取り出すシステムである。従来、この燃料電池の発電効率の向上を目的として、燃料電池内の化学反応解析や流動解析等を行うプログラムの開発がなされている。
【０００３】
例えば、特開２００５−４４６０２号公報においては、燃料電池の物理的及び化学的特性の解析方法並びにプログラムが提案されている（特許文献１）。この特許文献１に記載のプログラムは、はじめに初期値の設定を行い、その後、酸素濃度、電流密度、水素濃度、温度分布、化学反応をそれぞれ解析し、これらのデータをもとに燃料電池全体の電位が釣り合うか否かの判別をし、さらに電流値が実験値と一致するか否かの判別を行うものであり、電気化学的拘束条件下の物質や熱の流れを容易に、かつ精度よく解析することができると記載されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００５−４４６０２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、特許文献１に記載された発明を含む、これまでの燃料電池の数値解析プログラムは以下のような問題点を有している。すなわち、
１．従来の数値解析プログラムは、発電セルの平面内における所定の物理量の分布を解析するものである。よって、それぞれ性質の異なる多層領域間の物理量の移動を正確に求めているわけではなく、多層領域間を移動する物理量の計算や当該物理量の計算に基づいた発電量の計算等の解析結果がその精度において不十分であった。
２．燃料電池などの多層領域間における物理量は、各層の境界における化学反応等により各種の物理量の変化を伴い、解析結果は前記境界における境界条件に強い影響を受ける。そのため、与えた境界条件や初期値によっては収束性が悪くなったり、解が得られない場合があった。
３．燃料電池の解析において、多孔質材と電解質膜との間にある触媒層は、他の層に比べて非常に薄く、解析精度は触媒層近傍における格子解像度に依存していた。そのため、格子を十分に小さくする必要があったが、その分、計算負荷が大きくなっていた。
【０００６】
本発明は、このような問題点を解決するためになされたものであって、解析精度の向上、解析速度の向上、解析負荷の軽減ならびに収束性の安定を図ることのできる輸送方程式解析プログラムおよび輸送方程式解析方法を提供することを目的としている。
【課題を解決するための手段】
【０００７】
本発明に係る輸送方程式解析方法を用いた輸送方程式解析プログラムは、それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析プログラムであって、前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出部と、前記境界流束生成量算出部により算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第二領域内解析手段は、前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出部と、前記境界変化量算出部により算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出部と、前記第二領域側境界物理量算出部により算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析部と、前記第二領域内物理量解析部により求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出部として機能させる。
【０００８】
すなわち、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法は、第一領域と第二領域との境界において、ディレクレ条件に基づく境界条件とノイマン条件に基づく境界条件の２つの境界条件式を用いるとともに、その２つの境界条件式を所定の順序で解くことによって、汎用的な輸送方程式解析プログラムおよび輸送方程式解析方法を構築したものである。
【０００９】
また、本発明に係る輸送方程式解析方法を用いた輸送方程式解析プログラムは、それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析プログラムであって、前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出部と、前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出部と、前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第三領域内解析手段は、前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出部と、前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出部と、前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出部と、前記第三領域側境界物理量算出部により算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析部と、前記第三領域内物理量解析部により求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出部として機能させる。
【００１０】
すなわち、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法は、第二領域が薄膜近似と捉えられる場合であって、前記第一領域側境界および前記第三領域側境界に、ディレクレ条件に基づく境界条件とノイマン条件に基づく境界条件の２つの境界条件式を用いるとともに、それらの合計４つの境界条件式を所定の順序で解くことによって、薄膜近似を用いた汎用的な輸送方程式解析プログラムおよび輸送方程式解析方法を構築したものである。
【発明の効果】
【００１１】
本発明によれば、解析精度の向上、解析速度の向上、解析負荷の軽減ならびに収束性の安定を図ることができる。
【図面の簡単な説明】
【００１２】
【図１】本発明に係る輸送方程式解析プログラムを備えたコンピュータの第一実施形態を示すブロック図である。
【図２】本第一実施形態における輸送方程式解析プログラムにより解析を行う第一領域および第二領域を示す概略図である。
【図３】本第一実施形態の輸送方程式解析プログラムおよび輸送方程式解析方法における処理の手順を示すフローチャート図である。
【図４】本第一実施形態における第一領域内解析手段の処理の手順を示すフローチャート図である。
【図５】本第一実施形態における第二領域内解析手段の処理の手順を示すフローチャート図である。
【図６】本第一実施形態における各構成の機能を示す機能ブロック図である。
【図７】本発明に係る輸送方程式解析プログラムを備えたコンピュータの第二実施形態を示すブロック図である。
【図８】本第二実施形態における輸送方程式解析プログラムにより解析を行う第一領域ないし第三領域を示す概略図である。
【図９】本第二実施形態の輸送方程式解析プログラムおよび輸送方程式解析方法における処理の手順を示すフローチャート図である。
【図１０】本第二実施形態における第一領域内解析手段の処理の手順を示すフローチャート図である。
【図１１】本第二実施形態における第三領域内解析手段の処理の手順を示すフローチャート図である。
【図１２】本第二実施形態における各構成の機能を示す機能ブロック図である。
【図１３】高分子形燃料電池の概略構成を示す概略構成図でる。
【図１４】本実施例１における電解膜、触媒層および多孔質層での液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λの輸送方程式および境界条件を示す解析場概略図である。
【図１５】本実施例２における電解膜、触媒層および多孔質層での液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λの輸送方程式および境界条件を示す解析場概略図である。
【発明を実施するための形態】
【００１３】
以下、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法の第一実施形態について図面を用いて説明する。図１は、本第一実施形態における輸送方程式解析プログラム１ａを備えたコンピュータの構成を示すブロック図である。
【００１４】
また、図２は、本第一実施形態で解析を行う第一領域および第二領域を示す概略図である。前記第一領域および前記第二領域は、境界を介して連結されており、それぞれ性質の異なる構造および／または物質からなる。また、前記各領域の間を移動する所定の物理量は、異なる輸送方程式により表される。例えば、第一領域における輸送方程式は、以下の式（１）により表される。
・・・式（１）
また、第二領域における輸送方程式は、以下の式（２）により表される。
・・・式（２）
ここで、上記式（１）における「φ」は第一領域における物理量であり、上記式（２）における「ψ」は第二領域における物理量である。なお、本発明において、上記「φ」と上記「ψ」は同じ物理量を表すものであり、上記式（１）および上記式（２）における物理量を同じ符号を用いて記載することも可能であるが、以下に記載される各式における物理的意味を明確にするため、第一領域における物理量と第二領域における物理量を別表記とした。
【００１５】
また、本発明における物理量は、例えば、水、水蒸気および熱等を示すものであり、各領域において輸送方程式により表されるもの全ての物理量である。
【００１６】
本第一実施形態におけるコンピュータ１Ａは、図１に示すように、主として、入力手段２、表示手段３、記憶手段４および演算処理手段５から構成されている。以下、各構成について詳細に説明する。
【００１７】
入力手段２は、テキストや数値、記号等を入力する操作キー、操作マウス等からなる。本第一実施形態においては、輸送方程式解析プログラム１ａを実行する際に用いられる初期条件や設定条件等の数値の入力操作等に用いることができるようになっている。
【００１８】
表示手段３は、画像やテキストデータを表示する液晶ディスプレーやＣＲＴディスプレー、タッチパネル等からなり、輸送方程式解析プログラム１ａにおけるユーザーインターフェースとして、入力手段２により入力された内容や解析結果等を表示できるようになっている。
【００１９】
記憶手段４は、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、ハードディスク、フラッシュメモリ等によって構成されており、各種のデータを記憶するとともに、演算処理手段５が演算を行う際のワーキングエリアとして機能するものである。
【００２０】
本第一実施形態において、記憶手段４は、図１に示すように、主として、輸送方程式解析プログラム１ａを記憶するプログラム記憶部４１と、空間的に離散化された各格子点における物理量を記憶させておく物理量記憶部４２とを有する。なお、本第一実施形態における物理量記憶部４２には、各格子点における物理量とともに、各格子点における流束が記憶されている。
【００２１】
演算処理手段５は、ＣＰＵ（Central Processing Unit）等から構成されており、記憶手段４にインストールされた第一実施形態の輸送方程式解析プログラム１ａを実行させることにより、図１に示すように、物理量取得手段６、第一領域内解析手段７および第二領域内解析手段８としてコンピュータ１Ａを機能させるようになっている。以下、各構成手段についてより詳細に説明する。
【００２２】
物理量取得手段６は、記憶手段４の物理量記憶部４２に記憶された所定の時間ｔの各格子点における物理量を取得するものである。また、本第一実施形態における物理量取得手段６は、各格子点における物理量とともに、各格子点の流束も取得している。なお、流束は取得した物理量から計算して求めることも可能である。
【００２３】
第一領域内解析手段７は、物理量取得手段６により取得された物理量の時間ｔにより、所定時間Δｔ後における第一領域内の物理量の数値解を求めるものである。例えば、物理量取得手段６により取得された物理量が初期値である場合、物理量取得手段６により取得された物理量の時間ｔは、ｔ＝０であり、第一領域内解析手段７は、その時間ｔ＝０から所定時間Δｔ後の物理量を求めるものである。
【００２４】
本第一実施形態における第一領域内解析手段７は、図１に示すように、境界流束生成量算出部７１と、第一領域側境界流束算出部７２と、第一領域内物理量解析部７３と、第一領域側境界物理量算出部７４とを有している。
【００２５】
境界流束生成量算出部７１は、以下の式（３）の流束境界条件式により定義された境界流束生成量γに係る関数γ（φ，ψ）を用いて前記境界流束生成量γの値を算出するものである。
・・・（３）
【００２６】
よって、境界流束生成量算出部７１では、前記境界流束生成量γの関数γ（φ，ψ）に、前記物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界流束生成量γを算出している。
【００２７】
なお、この境界流束生成量γは、境界における化学反応等により物理量が変化したことにより生じる流束の差を示す値であり、算出された境界流束生成量γの符号が正の場合には単位時間当たりの生成量を、符号が負の場合には単位時間当たりの消滅量を表す。
【００２８】
また、上記式（３）で示される流束境界条件式は、第一領域と第二領域との間のノイマン条件による境界条件を表す式であり、前記第一領域側から前記境界へ流入出する前記物理量の流束ｆ（φ）と前記第二領域側から前記境界へ流入出する前記物理量の流束ｇ（ψ）との連続性を表したものである。
【００２９】
ここで、「α」は第一領域における係数であり、「β」は第二領域における係数である。例えば、物理量が液水であって第一領域が多孔質層である場合、液水は多孔質層内の孔の中しか通過することができない。よって、第一領域における流束は、流束ｎｆ（φ）と第一領域の多孔率（孔の体積／全体の体積）εとの積で表される。よって、この場合はα＝εとなる。
【００３０】
第一領域側境界流束算出部７２は、前記境界流束生成量算出部７１により算出された境界流束生成量γ、および前記物理量取得手段により取得された第二領域における境界近傍の流束ｇ＋（ψｋ）を、流束境界条件式に代入して所定時間Δｔ経過後の第一領域における境界近傍の流束ｆ−（φｋ＋１）を算出するものである。よって、第一領域側境界流束算出部７２における流束境界条件式は、所定時間Δｔ経過後の流束ｆ−（φｋ＋１）を算出するため、以下の式（４）に示すように時間方向に離散化されている。
・・・式（４）
【００３１】
ここで、「ｋ」は時間を表すものであって物理量取得手段６により取得された物理量の有する時間ｔに対応するものである。また、「ｋ＋１」は、時間ｔから所定時間Δｔ経過後の時間を表すものである。
【００３２】
また、「＋」は第二領域における境界近傍の値であることを示しており、具体的には、境界から第二領域側にある一つ目の格子点における値を示している。また、「−」は第一領域における境界近傍の値であることを示しており、具体的には、境界から第一領域側にある一つ目の格子点における値を示している。
【００３３】
よって、式（４）は、第二領域から境界へと流入する時間ｋにおける流束ｇ＋（ψｋ）が、境界において境界流束生成量γ分だけ増減され、所定時間Δｔ経過後の時間ｋ＋１には、第一領域における境界近傍へとｆ−（φｋ＋１）として流出していることを示している。
【００３４】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２により算出された境界近傍の流束αｎｆ−（φｋ＋１）を第一領域の輸送方程式に代入して所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求めるものである。よって、第一領域内物理量解析部７３における第一領域の輸送方程式は、所定時間Δｔ経過後の第一領域内における物理量の数値解を求めるために空間的に離散化すると、上記式（１）を以下の式（５）として表される。
・・・（５）
ここで、ｉは、第一領域内の格子点の位置を表すものであり、本第一実施形態においては、境界から２番目の格子点から、ｉ＝１，２，３，・・・と番号が割り当てられる。
【００３５】
本第一実施形態における第一領域内物理量解析部７３は、第一領域側境界流束算出部７２により算出した第一領域における境界近傍の流束ｆ−（φｋ＋１）を上記式（５）に代入して導かれる連立方程式を解くことで、第一領域内の各格子点における所定時間Δｔ経過後の物理量φｉｋ＋１を算出する。なお、上記式（５）による連立方程式を解く方法は特に限定されるものではなく、陰解法や陽解法等から適宜選択されるものである。
【００３６】
第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の物理量ａφ−ｋ＋１を算出するものである。
【００３７】
次に、第二領域内解析手段８について説明する。第二領域内解析手段８は、所定時間Δｔ経過後における第二領域内の物理量の数値解を求めるものである。本第一実施形態における第二領域内解析手段８は、図１に示すように、境界変化量算出部８１と、第二領域側境界物理量算出部８２と、第二領域内物理量解析部８３と、第二領域側境界流束算出部８４とを有している。
【００３８】
境界変化量算出部８１は、以下の式（６）の物理量境界条件式により定義された境界変化量ｃに係る関数ｃ（φ，ψ）を用いて前記境界変化量ｃの値を算出するものである。
・・・（６）
【００３９】
よって、境界変化量算出部８１では、前記境界変化量ｃの関数ｃ（φ，ψ）に、前記物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界変化量ｃを算出している。なお、この境界変化量ｃは、境界における流束の伝達率により生じる第一領域側の物理量と第二領域側の物理量との差を示す値である。
【００４０】
また、上記式（６）で示される物理量境界条件式は、第一領域と第二領域との間のディレクレ条件による境界条件を表す式であり、境界における第一領域側の物理量φと第二領域側の物理量ψとの連続性を示している。
【００４１】
第二領域側境界物理量算出部８２と境界変化量算出部８１とにより算出された境界変化量ｃおよび第一領域における境界近傍の物理量φ−ｋを物理量境界条件式に代入して、所定時間Δｔ経過後の第二領域における境界近傍の物理量ψ＋ｋ＋１を算出するものである。よって、第二領域側境界物理量算出部８２における物理量境界条件式は、所定時間Δｔ経過後の物理量ｂ（ψｋ）ψ＋ｋ＋１を算出するため、以下の式（７）に示すように時間方向に離散化されている。
・・・（７）
【００４２】
ここで、「ａ」は第一領域における係数であり第一領域の物理量φの関数として与えられている。また、「ｂ」は第二領域における係数であり第二領域の物理量ψの関数として与えられる。
【００４３】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２により算出された境界近傍の物理量ｂ（ψｋ）ψ＋ｋ＋１を第二領域の輸送方程式に代入して所定時間Δｔ経過後の第二領域内における物理量ψｊｋ＋１の数値解を求めるものである。よって、第二領域内物理量解析部８３における第二領域の輸送方程式は、所定時間Δｔ経過後の第二領域内における物理量ψｊｋ＋１の数値解を求めるために空間的に離散化すると、上記式（２）を以下の式（８）として表される。
・・・（８）
ここで、ｊは、第二領域内の格子点の位置を表すものであり、本第一実施形態においては、境界から２番目の格子点から、ｊ＝１，２，３，・・・と番号が割り当てられる。
【００４４】
本第一実施形態における第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２により算出した第二領域における境界近傍の物理量ψ＋ｋ＋１を上記式（８）に代入して導かれる連立方程式を解くことで、第二領域内の各格子点における所定時間Δｔ経過後の物理量ψｊｋ＋１を算出している。なお、上記式（８）による連立方程式を解く方法は特に限定されるものではなく、陰解法や陽解法等から適宜選択されるものである。
【００４５】
第二領域側境界流束算出部８４は、第二領域内物理量解析部８３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の流束ｇ＋（ψｋ＋１）を算出するものである。
【００４６】
次に、本第一実施形態の輸送方程式解析プログラム１ａにおける各構成の作用および輸送方程式解析方法について、図３ないし図５に示すフローチャート図および図６示す各構成の機能を示す機能ブロック図を用いて説明する。
【００４７】
図３および図６に示すように、本第一実施形態の輸送方程式解析プログラム１ａおよび輸送方程式解析方法では、まず、物理量取得手段６が、物理量記憶部４２に記憶された各格子点における物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を取得する（ステップＳ１）。続いて、サブルーチン化された第一領域内解析手段７が、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求める（ステップＳ２）。続いて、サブルーチン化された第二領域内解析手段８は、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ後の第二領域内における物理量ψｊｋ＋１の数値解を求める（ステップＳ３）。そして、時間ｔが所望する時間Ｔに達したか否かが判別される（ステップＳ４）。
【００４８】
なお、本第一実施形態では、図３に示すように、第一領域内解析手段７の処理後に第二領域内解析手段８の処理を行っているが、どちらを先に処理してもよく、並列に処理してもよい。また、本第一実施形態では、第一領域内解析手段７および第二領域内解析手段８をサブルーチン化して解析処理を実行しているが、サブルーチン化せずに物理量取得手段６と一連で処理するようにしてもよい。
【００４９】
以下、第一領域内解析手段７および第二領域内解析手段８について、より詳細に説明する。
【００５０】
図４に示すように、ステップＳ２における第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（３）により表される流束境界条件式において定義される境界流束生成量γに係る関数γ（φ，ψ）に、物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界流束生成量γの算出を行う（ステップＳ２１）。算出された境界流束生成量γは、図６に示すように、次のステップの第一領域側境界流束算出部７２に送られる。
【００５１】
第一領域側境界流束算出部７２は、時間方向に離散化された上記式（４）により表される流束境界条件式に、境界流束生成量算出部７１により算出された境界流束生成量γ、および物理量取得手段６により取得された第二領域における境界近傍の流束ｇ＋（ψｋ）を代入して、所定時間Δｔ経過後の第一領域における境界近傍の流束ｆ−（φｋ＋１）を算出する（ステップＳ２２）。算出された第一領域における境界近傍の流束ｆ−（φｋ＋１）は、図６に示すように、次のステップの第一領域内物理量解析部７３に送られるとともに、物理量記憶部４２に送られる。
【００５２】
物理量記憶部４２に送られた所定時間Δｔ経過後の第一領域における境界近傍の流束ｆ−（φｋ＋１）は、第一領域における境界近傍の格子点における流束値として置き換えられる。このように、物理量記憶部４２内の物理量を所定時間Δｔ経過後の物理量に置き換えることにより、時間発展計算を行うことができる。
【００５３】
次に、第一領域内物理量解析部７３は、時間方向に離散化された上記式（５）により表される第一領域の輸送方程式に、第一領域側境界流束算出部７２により算出された境界近傍の流束ｆ−（φｋ＋１）を代入して導き出される連立方程式を解くことで、第一領域内の各格子点における所定時間Δｔ経過後の物理量φｉｋ＋１の算出を行う（ステップＳ２３）。算出された第一領域における物理量φｉｋ＋１は、図６に示すように、次のステップの第一領域側境界物理量算出部７４に送られるとともに、物理量記憶部４２に送られる。
【００５４】
第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の物理量φ−ｋ＋１を算出する（ステップＳ２４）。そして、算出された第一領域における境界近傍の物理量φ−ｋ＋１は、図６に示すように、物理量記憶部４２に送られる。以上で、第一領域内解析手段７によるステップＳ２は終了し、次の第二領域内解析手段８によるステップＳ３へと進む。
【００５５】
図５に示すように、ステップＳ３における第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（６）の物理量境界条件式により定義された境界変化量ｃに係る関数ｃ（φ，ψ）に、物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界変化量ｃの算出を行う（ステップＳ３１）。算出された境界変化量ｃは、図６に示すように、次のステップの第二領域側物理量算出部８２に送られる。
【００５６】
第二領域側境界物理量算出部８２は、時間方向に離散化された上記式（７）により表される物理量境界条件式に、境界変化量算出部８１により算出された境界変化量ｃ、および物理量取得手段６により取得された第一領域における境界近傍の物理量φ−ｋを代入して、所定時間Δｔ経過後の第二領域における境界近傍の物理量ψ＋ｋ＋１を算出する（ステップＳ３２）。算出された第二領域における境界近傍の物理量ψ＋ｋ＋１は、図６に示すように、次のステップの第二領域内物理量解析部８３に送られるとともに、物理量記憶部４２に送られる。
【００５７】
次に、第二領域内物理量解析部８３は、時間方向に離散化された上記式（８）により表される第二領域の輸送方程式に、第二領域側境界物理量算出部８２により算出された境界近傍の物理量ψ＋ｋ＋１を代入して導き出された連立方程式を解くことによって第二領域内の各格子点における所定時間Δｔ経過後の物理量ψｊｋ＋１の算出を行う（ステップＳ３３）。算出された第二領域における物理量ψｊｋ＋１は、図６に示すように、次のステップの第二領域側境界流束算出部８４に送られるとともに、物理量記憶部４２に送られる。
【００５８】
第一領域側境界流束算出部８４は、第二領域内物理量解析部８３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の流束ｇ＋（ψｋ＋１）を算出する（ステップＳ３４）。そして、算出された第二領域における境界近傍の流束ｇ＋（ψｋ＋１）は、図６に示すように、物理量記憶部４２に送られる。以上で、第二領域内解析手段８によるステップＳ３は終了し、次のステップＳ４へと進む。
【００５９】
ステップＳ４では、時間ｔが所望する時間Ｔに達したか否かを行う。時間ｔが所望する時間Ｔに達していない場合（ステップＳ４：ＮＯ）、時間ｔ＝ｔ＋ΔｔとしてステップＳ１へもどり、ステップ１からステップＳ３の処理を繰り返し行う。このように、繰り返し処理を行うことにより第一領域内および第二領域内の各格子点における物理量φ，ψの時間発展計算が行われる。そして、時間ｔが所望する時間Ｔに達した場合（ステップＳ４：ＹＥＳ）、処理を終了させる。
【００６０】
以上のような本第一実施形態の輸送方程式解析プログラム１ａおよび輸送方程式解析方法によれば、以下のような効果を得ることができる。
１．多層領域間を移動する物理量の数値解析を行うことができる。
２．流束境界条件式および物理量境界条件式の２つの境界条件を算出する式を用いることにより、適当な各領域間の境界における物理量を求めることができるため、各領域における輸送方程式の数値解の精度を向上させるとともに、収束性も向上する。
３．第一領域側と第二領域側との解析を並列処理することができるため、処理時間を向上させることができる。
４．２つの輸送方程式の連成問題であれば物理量の種類については特に限定がなく、応用範囲が広い。
【００６１】
次に、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法の第二実施形態について説明する。なお、本第二実施形態のうち、前述した第一実施形態の構成と同一若しくは相当する構成については同一の符号を付して再度の説明を省略する。
【００６２】
図７は、本第二実施形態における輸送方程式解析プログラム１ｂを備えたコンピュータの構成を示すブロック図である。
【００６３】
また、図８は、本第二実施形態で解析を行う第一領域、第二領域および第三領域を示す概略図である。図８に示すように、本第二実施形態における第一領域、第二領域および第三領域は、第一領域と第二領域とがこれらの間の第一領域側境界を介して連結されており、第二領域と第三領域とがこれらの間の第三領域側境界を介して連結されている。
【００６４】
また、第一領域および第三領域は、それぞれ性質の異なる構造および／または物質からなり、前記各領域の間を移動する所定の物理量は、異なる輸送方程式により表されるものである。具体的には、第一領域における輸送方程式は、上記式（１）により表されるととともに、第三領域における輸送方程式は、上記式（２）により表される。
【００６５】
さらに、本第二実施形態における第二領域は薄膜状に形成されている。そのため、第二領域内における物理量および流束は、空間的にほぼ一定であると見なすことができる。よって、本第二実施形態では、第二領域内の物理量およびその物理量の流束を、それぞれθおよびｈ（θ）と近似した。本発明では、以上のような第二領域における物理量等の近似を薄膜近似と定めている。
【００６６】
以下、本第二実施形態の各構成について、より詳細に説明する。
【００６７】
本第二実施形態におけるコンピュータ１Ｂは、図７に示すように、主として、入力手段２、表示手段３、記憶手段４および演算処理手段５から構成されている。
【００６８】
このうち本第二実施形態における演算処理手段５は、記憶手段４にインストールされた第二実施形態の輸送方程式解析プログラム１ｂを実行させることにより、図７に示すように、物理量取得手段６、第一領域内解析手段９および第三領域内解析手段１０としてコンピュータ１Ｂを機能させるようになっている。
【００６９】
本第二実施形態における第一領域内解析手段９は、上述した第一実施形態における第一領域内解析手段７と同様に、所定時間Δｔ経過後における第一領域内の物理量の数値解を求めるものである。
【００７０】
また、本第二実施形態における第一領域内解析手段９は、図７に示すように、第一領域側境界流束生成量算出部９１と、第三領域側境界流束生成量算出部９２と、第一領域側境界流束算出部９３と、第一領域内物理量解析部９４と、第一領域側境界物理量算出部９５とを有している。つまり、上述した第一実施形態では、境界が１つしかないため境界の境界流束生成量の算出は、境界流束生成量算出部７１による１つの処理で行うことができるのに対し、本第二実施形態では、境界が、第一領域側境界と第三領域側境界の２つあるため、各境界における流束の算出についても第一領域側境界流束生成量算出部９１と第三領域側境界流束生成量算出部９２とによる２つの処理を要する。
【００７１】
第一領域側境界流束生成量算出部９１は、以下の式（９）の流束境界条件式により定義された境界流束生成量γ−に係る関数γ−（φ，θ）を算出するものであり、前記関数γ−（φ，θ）に取得された物理量φおよび物理量θを代入して第一領域側境界における物理量の境界流束生成量γ−を算出している。
・・・（９）
【００７２】
また、上記式（９）で示される流束境界条件式は、第一領域と第二領域との間のノイマン条件による境界条件を表す式であり、第一領域側から第一領域側境界へ流入出する物理量の流束ｆ（φ）と第二領域側から第一領域側境界へ流入出する前記物理量の流束ｈ（θ）との連続性を示している。ここで、「−」は、第一領域側境界の値であることを示している。
【００７３】
次に、第三領域側境界流束生成量算出部９２は、以下の式（１０）の流束境界条件式により定義された境界流束生成量γ＋に係る関数γ＋（θ，ψ）を算出するものであり、前記関数γ＋（θ，ψ）に物理量取得手段６により取得された物理量ψおよびθを代入して第三領域側境界における物理量の境界流束生成量γ＋を算出している。
・・・（１０）
【００７４】
また、上記式（１０）で示される流束境界条件式は、第三領域と第二領域との間のノイマン条件による境界条件を表す式であり、第三領域側から第三領域側境界へ流入出する物理量の流束ｇ（ψ）と第二領域側から第三領域側境界へ流入出する前記物理量の流束ｈ（θ）との連続性を示している。ここで、「＋」は、第三領域側境界の値であることを示している。
【００７５】
第一領域側境界流束算出部９３は、第一領域側境界流束生成量算出部９１により算出された境界流束生成量γ−、第三領域側境界流束生成量算出部９２により算出された境界流束生成量γ＋、および前記物理量取得手段６により取得された第三領域側境界の流束ｇ＋（ψｋ）を、前記流束境界条件式に代入して所定時間Δｔ経過後の第一領域側境界の流束ｆ−（φｋ＋１）を算出するものである。そのため、第一領域側境界流束算出部９３では、時間方向に離散化された、以下の式（１１）および式（１２）からなる流束境界条件式が用いられる。
・・・（１１）
・・・（１２）
【００７６】
具体的には、まず、上記式（１１）に第三領域側境界流束生成量算出部９２により算出された境界流束生成量γ＋と、第三領域境界の流束ｇ＋（ψｋ）を代入して、所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を算出する。そして、上記式（１２）に、第一領域側境界流束生成量算出部９１により算出された境界流束生成量γ−と式（１１）により算出された所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を代入して、所定時間Δｔ経過後の第一領域側境界の流束ｆ−（φｋ＋１）を算出する。
【００７７】
本第二実施形態における第一領域内物理量解析部９４は、第一実施形態における第一領域内物理量解析部７３と同様の機能を有しており、第一領域側境界流束算出部９３により算出された第一領域側境界の流束ｆ−（φｋ＋１）を、上記式（５）で表される第一領域の輸送方程式に代入して所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求めるものである。
【００７８】
また、本第二実施形態における第一領域側境界物理量算出部９５は、第一実施形態における第一領域側境界物理量算出部７４と同様の機能を有しており、第一領域内物理量解析部９４により求められた数値解から所定時間Δｔ経過後の第一領域側境界の物理量φ−ｋ＋１を算出するものである。
【００７９】
次に、第三領域内解析手段１０について説明する。本第二実施形態における第三領域内解析手段１０は、第一実施形態における第二領域内解析手段８と同様に、所定時間Δｔ経過後における第三領域内の物理量の数値解を求めるものである。本第二実施形態における第二領域内解析手段１０は、図７に示すように、第三領域側境界変化量算出部１０１と、第一領域側境界変化量算出部１０２と、第三領域側境界物理量算出部１０３と、第三領域内物理量解析部１０４と、第三領域側境界流束算出部１０５とを有している。
【００８０】
第三領域側境界変化量算出部１０１は、以下の式（１３）の流束境界条件式により定義された境界変化量ｃ＋に係る関数ｃ＋（θ，ψ）を算出するものであり、前記関数ｃ＋（θ，ψ）に前記物理量取得手段６により取得された物理量ψおよび物理量θを代入して第三領域側境界における物理量の境界変化量ｃ＋を算出している。
・・・（１３）
【００８１】
また、上記式（１３）で示される流束境界条件式は、第三領域と第二領域との間のディレクレ条件による境界条件を表す式であり、第三領域側境界における第三領域側の物理量ψと第二領域側の物理量θとの連続性を示している。
【００８２】
第一領域側境界変化量算出部１０２は、以下の式（１４）の流束境界条件式により定義された境界変化量ｃ−に係る関数ｃ−（φ，θ）を算出するものであり、前記関数ｃ−（φ，θ）に前記物理量取得手段６により取得された物理量φおよび物理量θを代入して第一領域側境界における物理量の境界変化量ｃ−を算出している。
・・・（１４）
【００８３】
また、上記式（１４）で示される流束境界条件式は、第一領域と第二領域との間のディレクレ条件による境界条件を表す式であり、第一領域側境界における第一領域側の物理量φと第二領域側の物理量θとの連続性を示している。
【００８４】
第三領域側境界物理量出部１０３は、第三領域側境界変化量算出部１０１により算出された境界変化量ｃ＋、第一領域側境界変化量算出部１０２により算出された境界変化量ｃ−、および前記物理量取得手段６により取得された第一領域境界の物理量φ−を、物理量境界条件式に代入して所定時間Δｔ経過後の第三領域側境界の物理量ψ＋ｋ＋１を算出するものである。そのため、第三領域側境界物理量算出部１０３では、時間方向に離散化された、以下の式（１５）および式（１６）からなる物理量境界条件式が用いられる。
・・・（１５）
・・・（１６）
【００８５】
具体的には、まず、上記式（１５）に第一領域側境界変化量算出部１０２により算出された境界変化量ｃ−と、第一領域境界の物理量φ−を代入して、所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を算出する。そして、上記式（１６）に、第三領域側境界変化量算出部１０１により算出された境界変化量ｃ＋と、式（１５）により算出された所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を代入して、所定時間Δｔ経過後の第三領域側境界の物理量ψ−ｋ＋１を算出する。
【００８６】
本第二実施形態における第三領域内物理量解析部１０４は、第一実施形態における第二領域内物理量解析部８３と同様の機能を有しており、第三領域側境界物理量算出部１０３により算出された第三領域側境界の物理量ψ−ｋ＋１を、上記式（８）で表される第三領域の輸送方程式に代入して所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求めるものである。
【００８７】
また、本第二実施形態における第三領域側境界流束算出部１０５は、第一実施形態における第二領域側境界流束算出部８４と同様の機能を有しており、第三領域内物理量解析部１０４により求められた数値解から所定時間Δｔ経過後の第三領域側境界の物理量の流束ｇ＋（ψｋ＋１）を算出するものである。
【００８８】
次に、本第二実施形態の輸送方程式解析プログラム１ｂにおける各構成の作用および輸送方程式解析方法について、図９ないし図１１に示すフローチャート図および図１２に示す各構成の機能を示す機能ブロック図を用いて説明する。
【００８９】
図９および図１２に示すように、本第二実施形態の輸送方程式解析プログラム１ｂおよび輸送方程式解析方法では、まず、物理量取得手段６が、物理量記憶部４２に記憶された各格子点における物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を取得する（ステップＳＳ１）。続いて、サブルーチン化された第一領域内解析手段９が、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求める（ステップＳＳ２）。続いて、サブルーチン化された第三領域内解析手段１０は、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ経過後の第三領域内における物理量ψｊｋ＋１の数値解を求める（ステップＳＳ３）。そして、時間ｔが所望する時間Ｔに達したか否かが判別される（ステップＳＳ４）。
【００９０】
なお、第一領域内解析手段９の処理と第三領域内解析手段１０の処理とはどちらを先に処理してもよく、並列に処理してもよい。また、本第二実施形態では、第一領域内解析手段９および第三領域内解析手段１０をサブルーチンとせずに物理量取得手段６と一連で処理するようにしてもよい。
【００９１】
図１０に示すように、ステップＳＳ２における第一領域内解析手段９では、まず、第一領域側境界流束生成量算出部９１が、上記式（９）の第一領域側境界における流束境界条件式により定義された境界流束生成量γ−に係る関数γ−（φ，θ）に、物理量取得手段６により取得された物理量φおよび物理量θを代入して、第一領域側境界における物理量の境界流束生成量γ−の算出を行う（ステップＳＳ２１）。算出された境界流束生成量γ−は、図１２に示すように、第一領域側境界流束算出部９３に送られる。
【００９２】
また、第三領域側境界流束生成量算出部９２が、上記式（１０）により表される第三領域側境界における流束境界条件式により定義され境界流束生成量γ＋に係る関数γ＋（θ，ψ）に、物理量取得手段６で取得された物理量ψおよびθを代入して、第一領域側境界における物理量の境界流束生成量γ＋の算出を行う（ステップＳＳ２２）。算出された境界流束生成量γ＋は、図１２に示すように、第一領域側境界流束算出部９３に送られる。
【００９３】
第一領域側境界流束算出部９３は、時間方向に離散化された上記式（１１）により表される第三領域側境界における流束境界条件式に、第三領域側境界流束生成量算出部９２により算出された境界流束生成量γ＋、および物理量取得手段６により取得された第三領域側境界の流束ｇ＋（ψｋ）を代入して、所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を算出するとともに、上記式（１２）に、第一領域側境界流束生成量算出部９１により算出された境界流束生成量γ−と式（１１）により算出された所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を代入して、所定時間Δｔ経過後の第一領域側境界の流束ｆ−（φｋ＋１）を算出する（ステップＳＳ２３）。算出された第一領域側境界の流束ｆ−（φｋ＋１）は、図１２に示すように、次のステップの第一領域内物理量解析部９４に送られるとともに、物理量記憶部４２に送られる。
【００９４】
次に、第一領域内物理量解析部９４は、時間方向に離散化された上記式（５）により表される第一領域の輸送方程式に、第一領域側境界流束算出部９３により算出された第一領域側境界の流束ｆ−（φｋ＋１）を代入して導き出される連立方程式を解くことにより、第一領域内の各格子点における所定時間Δｔ経過後の物理量φｉｋ＋１の算出を行う（ステップＳＳ２４）。算出された第一領域における物理量φｉｋ＋１は、図１２に示すように、次のステップの第一領域側境界物理量算出部９５に送られるとともに、物理量記憶部４２に送られる。
【００９５】
第一領域側境界物理量算出部９５は、第一領域内物理量解析部９４により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の物理量φ−ｋ＋１を算出する（ステップＳＳ２５）。そして、算出された第一領域における境界近傍の物理量φ−ｋ＋１は、図１２に示すように、物理量記憶部４２に送られる。
【００９６】
図１１に示すように、ステップＳＳ３における第二領域内解析手段１０では、まず、第三領域側境界変化量算出部１０１が、上記式（１３）の第三領域側境界における物理量境界条件式により定義された境界変化量ｃ＋（θ，ψ）に係る関数に、物理量取得手段６により取得された物理量ψおよび物理量θを代入して、第三領域側境界における物理量の境界変化量ｃ＋の算出を行う（ステップＳＳ３１）。算出された境界変化量ｃ＋は、図１１に示すように、第三領域側物理量算出部１０３に送られる。
【００９７】
次に、第一領域側境界変化量算出部１０２が、上記式（１４）の第一領域側境界における物理量境界条件式により定義された境界変化量ｃ−（φ，θ）に係る関数に、物理量取得手段６により取得された物理量φおよび物理量θを代入して、第一領域側境界における物理量の境界境界変化量ｃ−の算出を行う（ステップＳＳ３２）。算出された境界変化量ｃ−は、図１１に示すように、第三領域側物理量算出部１０３に送られる。
【００９８】
第三領域側境界物理量算出部１０３は、時間方向に離散化された上記式（１５）により表される第一領域側境界における流束境界条件式に、第一領域側境界変化量算出部１０２により算出された境界変化量ｃ−、および物理量取得手段６により取得された第一領域側境界の物理量φｋを代入して、所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を算出するとともに、上記式（１６）に、第三領域側境界変化量算出部１０１により算出された境界変化量ｃ＋と、式（１５）により算出された所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を代入して、所定時間Δｔ経過後の第三領域側境界の物理量ψ＋ｋ＋１を算出する（ステップＳＳ２３）。算出された第三領域側境界の物理量ψ＋ｋ＋１は、図１２に示すように、次のステップの第三領域内物理量解析部１０４に送られるとともに、物理量記憶部４２に送られる。
【００９９】
次に、第三領域内物理量解析部１０４は、時間方向に離散化された上記式（８）により表される第三領域の輸送方程式に、第三領域側境界物理量算出部１０３により算出された第三領域側境界における物理量ψ＋ｋ＋１を代入して導き出される連立方程式を解くことにより、第三領域内の各格子点における所定時間Δｔ経過後の物理量ψｊｋ＋１の算出を行う（ステップＳＳ３４）。算出された第三領域側境界における物理量ψｊｋ＋１は、図１２に示すように、次のステップの第三領域側境界流束算出部１０５に送られるとともに、物理量記憶部４２に送られる。
【０１００】
第三領域側境界流束算出部１０５は、第三領域内物理量解析部１０４により求められた数値解から所定時間Δｔ経過後の第三領域側境界における流束ｇ＋（ψｋ＋１）を算出する（ステップＳＳ３５）。そして、算出された第三領域側境界の流束ｇ＋（ψｋ＋１）は、図１２に示すように、物理量記憶部４２に送られる。
【０１０１】
ステップＳＳ４では、時間ｔが所望する時間Ｔに達したか否かを行う。時間ｔが所望する時間Ｔに達していない場合（ステッＳプＳ４：ＮＯ）、時間ｔ＝ｔ＋Δｔとし、ステップＳＳ１からステップＳＳ３の処理を繰り返し行う。そして、時間ｔが所望する時間Ｔに達した場合（ステップＳＳ４：ＹＥＳ）、処理を終了させる。
【０１０２】
以上のような本第二実施形態の輸送方程式解析プログラム１ｂおよび輸送方程式解析方法によれば、薄膜近似した第二領域にける格子解像度を高くする必要が無くなるため、計算負荷を小さくすることができるという効果が得られる。
【実施例１】
【０１０３】
『本発明を用いた高分子形燃料電池の解析』
本実施例１では、上述した第一実施形態の解析方法を用いて、高分子形燃料電池内における各種物理量の輸送方程式および各層の境界における境界条件式を例示しつつ、解析について説明する。
【０１０４】
高分子形燃料電池は、図１３に示すように、電解膜と、この電界膜を中心に陰極側および陽極側にそれぞれ触媒層と多孔質層とが多層状に連結されている。また、高分子形燃料電池では、水素や酸素などの作動ガス、生成物である水または水蒸気等が混合された状態で前記各層内を移動する、いわゆる複合流動場が構成されている。
【０１０５】
そこで、本実施例１では、各層における物理量を液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λに分けて、それぞれについて輸送方程式が設定されている。また、液水ｓおよび水蒸気質量分率Ｙｗは直接相変化が行われるとともに、水蒸気質量分率Ｙｗおよび電界質内の含水量λは直接相変化が行われるものとし、一方、液水ｓおよび電界質内の含水量λは直接相変化が行われないものとした。
【０１０６】
（１）各層における輸送方程式について
各層における輸送方程式について、図１４を参照しながら説明する。図１４に示すように、以下の式（１７）は、触媒層内における液水ｓの輸送方程式である。
・・・（１７）
また、以下の式（１８）は、多孔質層内における液水ｓの輸送方程式である。
・・・（１８）
【０１０７】
ここで、上記式（１７）および上記式（１８）内の各係数は、以下の式（１９）ないし式（２１）の関係を有している。
・・・（１９）
・・・（２０）
・・・（２１）
【０１０８】
また、εは多効率、ρｌは水の密度、ρｇはガスの密度、ｕ気体の速度、Ｄｃ拡散係数、κ透過率、νは動粘性計数、ｇは重力加速度、Ｍｌはモル質量、ｋｃは凝縮係数、ｋｅは蒸発係数、Ｘｗはモル分率、Ｒは気体定数、Ｔ温度、Ｐ圧力、Ｐｓａｔは飽和水蒸気量、σは表面張力、θｃは液水ｓの接触角およびμは粘性係数を示している。また、ｍ１は液水ｓと水蒸気質量分率Ｙｗの間の相変化量である。
【０１０９】
なお、電界膜は固体であり液水ｓは入り込めないため、図１４に示すように、電界膜内の液水ｓ＝０である。
【０１１０】
次に、水蒸気質量分率Ｙｗの輸送方程式について説明する。以下の式（２２）は、触媒層内における水蒸気質量分率Ｙｗの輸送方程式である。
・・・（２２）
また、以下の式（２３）は、多孔質層内における水蒸気質量分率Ｙｗの輸送方程式である。
・・・（２３）
【０１１１】
また、上記式（２２）および上記式（２３）内の各係数は、以下の式（２４）の関係を有している。
・・・（２４）
【０１１２】
ここで、Ｙｉは化学種ｉの質量分率、Ｄｉ，ｅｆｆは参照拡散係数、Ｄｉ０は拡散係数、Ｐ０は基準圧力、Ｔ０は基準温度、ωは化学反応による水または水蒸気の発生量、ｍ２は水蒸気質量分率Ｙｗと電界質内の含水量λと間の相変化量である。
【０１１３】
また、水蒸気質量分率Ｙｗは、液水ｓと同様に、電界膜内に入り込めないため、電界膜内では水蒸気質量分率Ｙｗ＝０である。
【０１１４】
次に、電界質内の含水量λの輸送方程式について説明する。以下の式（２５）は、電解膜における電界質内の含水量λの輸送方程式である。
・・・（２５）
また、以下の式（２６）は、触媒層内における電解質内の含水量λの輸送方程式である。
・・・（２６）
【０１１５】
また、上記式（２５）および上記式（２６）内の各係数は、以下の式（２７）および式（２８）の関係を有している。
・・・（２７）
・・・（２８）
【０１１６】
ここで、ρｍは電解膜の密度、ＥＷは電解膜の等価質量、Ｖｅｘは電解膜の膨張係数、Ｄｗｌは電解膜の拡散係数、Ｄｉは拡散係数、ｎｄは電気浸透抵抗係数、ｉは電流、Ｆはファラデー定数である。
【０１１７】
なお、電界膜は固体であり液水ｓは入り込めないため、電界膜内では液水ｓ＝０である。
【０１１８】
以上のように、高分子形燃料電池の問題においては、液水ｓは触媒層と多孔質層の二つの領域を移動する物理量であり、各領域にいて異なる輸送方程式により示すことができる。また、水蒸気質量分率Ｙｗは触媒層と多孔質層の二つの領域を移動する物理量であり、各領域において異なる輸送方程式により示すことができる。さらに、電界質内の含水量λは電解質と触媒層の二つの領域を移動する物理量であり、各領域において異なる輸送方程式により示すことができる。
【０１１９】
すなわち、高分子形燃料電池内の液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λは、それぞれ性質の異なる第一領域および第二領域が第一領域と第二領域との間の境界を介して連結されており、各々の領域を移動する物理量が第一領域および第二領域で異なる輸送方程式によって表されるものであるため、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法によって処理可能な問題である。
【０１２０】
（２）液水ｓの輸送方程式の解析について
まず、液水ｓの輸送方程式の解析について説明する。本実施例１では、第一領域を触媒層、第二領域を多孔質層とした。
【０１２１】
液水ｓの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式（２９）で表される。
・・・（２９）
また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式（３０）で表される。
・・・（３０）
ただし、第一領域側のＰｃ，ＣＬは、以下の式（３１）および式（３２）によって表されるものである。
・・・（３１）
・・・（３２）
【０１２２】
なお、Ｐｃはｓの三次の多項式であり、逆関数Ｐｃ−１を直接求めることが難しい。そこで、以下の方法でＰｃの線形化を行った。
【０１２３】
まず、ｓ＝ａにおいてテイラー展開を行うと、以下の式（３３）で表される一次多項式が得られる。
・・・（３３）
上記式（３３）をＰｃに適用すると、以下の式（３４）が得られる。
・・・（３４）
【０１２４】
以上を踏まえると、上記式（３１）は線形化され、以下の式（３５）として表される。
・・・（３５）
また、第二領域側のＰｃ，ＧＤＬにも適用すると、Ｐｃ，ＧＤＬは以下の式（３６）として表される。
・・・（３６）
ただし、上記式（３５）および上記式（３６）は以下の関係を有する。
【０１２５】
ここで、ａは前時刻のｓ（ａ＝ｓｋ−１）を用いることとすると、上記式（３０）は、以下の式（３７）として表される。
・・・（３７）
【０１２６】
第一領域側の輸送方程式である上記式（１７）を離散化すると、以下の式（３８）として表される。
・・・（３８）
また、第二領域側の輸送方程式である上記式（１８）を離散化すると、以下の式（３９）として表される。
・・・（３９）
なお、上記式（３８）および上記式（３９）は、時間発展計算をオイラー陽解法で解く場合の式である。
【０１２７】
次に、各境界条件についても離散化を行う。以下の式（４０）は、ノイマン条件における上記式（２９）を離散化したものである。
・・・（４０）
また、以下の式（４１）は、ディレクレ条件における上記式（３７）を離散化したものである。
・・・（４１）
【０１２８】
次に、以上の各式を用いて、第一領域および第二領域における液水ｓの輸送方程式の解析について、特に、第一領域内解析手段７および第二領域内解析手段８の解析手順に従って、以下に説明する。
【０１２９】
第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（２９）により境界流束生成量γを算出する。本件の場合は、図１４に示すように、境界における流束（Ｊｆｌｕｘ）は、化学反応等により変化しないと仮定しているため、境界流束生成量γ＝０である。
【０１３０】
次に、第一領域側境界流束算出部７２は、上記式（４０）を用いて第一領域の境界近傍の流束ｆ−（ｓＣＬｋ＋１）を算出する。なお、本件における流束は、各領域における多孔率εに影響されるため、α＝εＣＬであり、β＝εＧＤＬである。
【０１３１】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２より算出された流束ｆ−（ｓＣＬｋ＋１）を用いて上記式（３８）を解き、液水ｓＣＬｋ＋１を求める。
【０１３２】
そして、第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた液水ｓＣＬｋ＋１の値を用いて、第一領域の境界近傍における液水ｓの値ｓ−ｋ＋１を算出する。
【０１３３】
以上により、第一領域内の所定時間Δｔ経過後における液水ｓの値を全て求めることができる。
【０１３４】
次に、第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（３７）により境界における液水ｓの境界変化量ｃを算出する。本件では、境界におけるＰｃの値は化学反応等により変化しないとしているため、境界変化量ｃ＝０である。
【０１３５】
第二領域側境界物質量算出部８２は、上記式（４１）を用いて第二領域の境界近傍の液水ｓＣＬｋ＋１を算出する。なお、本件では、ａ＝１でありｂ＝１である。
【０１３６】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２より算出された液水ｓＣＬｋ＋１を用いて上記式（３９）を解き、液水ＧＤＬｋ＋１を求める。
【０１３７】
そして、第二領域側境界流束算出部は、第二領域内物理量解析部８３により求められた液水ｓＧＤＬｋ＋１の値を用いて、第二領域の境界近傍における液水ｓの流束ｇ＋（ｓＧＤＬｋ＋１）を算出する。
【０１３８】
以上により、第二領域内の所定時間Δｔ経過後における液水ｓの値を全て求めることができる。
【０１３９】
（３）水蒸気質量分率Ｙｗの輸送方程式の解析について
次に、水蒸気質量分率Ｙｗの輸送方程式の解析について説明する。本実施例１では、第一領域を触媒層、第二領域を多孔質層とした。
【０１４０】
水蒸気質量分率Ｙｗの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式（４２）で表される。
・・・（４２）
また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式（４３）で表される。
・・・（４３）
【０１４１】
次に、第一領域側の輸送方程式である上記式（２２）を離散化すると、以下の式（４４）として表される。
・・・（４４）
また、第二領域側の輸送方程式である上記式（２３）を離散化すると、以下の式（４５）として表される。
・・・（４５）
なお、上記式（４４）および上記式（４５）は、時間発展計算をオイラー陽解法で解く場合の式である。
【０１４２】
次に、各境界条件についても離散化を行う。以下の式（４６）は、ノイマン条件における上記式（４２）を離散化したものである。
・・・（４６）
また、以下の式（４７）は、ディレクレ条件における上記式（４３）を離散化したものである。
・・・（４７）
【０１４３】
次に、以上の各式を用いて、第一領域および第二領域における水蒸気質量分率Ｙｗの輸送方程式の解析について、特に、第一領域内解析手段７および第二領域内解析手段８の解析手順に従って、以下に説明する。
【０１４４】
第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（４２）により境界流束生成量γを算出する。本件の場合は、図１４に示すように、境界における流束（Ｊｆｌｕｘ）は、化学反応等により変化しないと仮定しているため、境界流束生成量γ＝０である。
【０１４５】
次に、第一領域側境界流束算出部７２は、上記式（４６）を用いて第一領域の境界近傍の流束ｆ−（ＹＣＬｋ＋１）を算出する。なお、本件における流束は、各領域における多孔率εに影響されるため、α＝εＣＬであり、β＝εＧＤＬである。
【０１４６】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２より算出された流束ｆ−（ｓＣＬｋ＋１）を用いて上記式（４４）を解き、水蒸気質量分率ＹＣＬｋ＋１を求める。
【０１４７】
そして、第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた水蒸気質量分率ＹＣＬｋ＋１の値を用いて、第一領域の境界近傍における水蒸気質量分率Ｙの値Ｙ−ｋ＋１を算出する。
【０１４８】
以上により、第一領域内の所定時間Δｔ経過後における水蒸気質量分率Ｙｗの値を全て求めることができる。
【０１４９】
次に、第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（４３）により境界における水蒸気質量分率Ｙｗの境界変化量ｃを算出する。本件では、境界におけるＹｗの値は化学反応等により変化しないと仮定しているため、境界変化量ｃ＝０である。
【０１５０】
第二領域側境界物質量算出部８２は、上記式（４７）を用いて第二領域の境界近傍の水蒸気質量分率ＹＣＬｋ＋１を算出する。なお、本件では、ａ＝１でありｂ＝１である。
【０１５１】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２より算出された水蒸気質量分率ＹＣＬｋ＋１を用いて上記式（４５）を解き、水蒸気質量分率ＹＧＤＬｋ＋１を求める。
【０１５２】
そして、第二領域側境界流束算出部は、第二領域内物理量解析部８３により求められた水蒸気質量分率ＹＧＤＬｋ＋１の値を用いて、第二領域の境界近傍における水蒸気質量分率Ｙｗの流束ｇ＋（ＹＧＤＬｋ＋１）を算出する。
【０１５３】
以上により、第二領域内の所定時間Δｔ経過後における水蒸気質量分率Ｙｗの値を全て求めることができる。
【０１５４】
（４）電界質内の含水量λの輸送方程式の解析について
次に、電界質内の含水量λの輸送方程式の解析について説明する。本実施例１では、第一領域を電解膜、第二領域を触媒層とした。
【０１５５】
電界質内の含水量λの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式（４８）で表される。
・・・（４８）
また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式（４９）で表される。
・・・（４９）
【０１５６】
次に、第一領域側の輸送方程式である上記式（２５）を離散化すると、以下の式（５０）として表される。
・・・（５０）
また、第二領域側の輸送方程式である上記式（２６）を離散化すると、以下の式（５１）として表される。
・・・（５１）
なお、上記式（５０）および上記式（５１）は、時間発展計算をオイラー陽解法で解く場合の式である。
【０１５７】
次に、各境界条件についても離散化を行う。以下の式（５２）は、ノイマン条件における上記式（４８）を離散化したものである。
・・・（５２）
また、以下の式（５３）は、ディレクレ条件における上記式（４９）を離散化したものである。
・・・（５３）
【０１５８】
次に、以上の各式を用いて、第一領域および第二領域における電界質内の含水量λの輸送方程式の解析について、特に、第一領域内解析手段７および第二領域内解析手段８の解析手順に従って、以下に説明する。
【０１５９】
第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（４８）により境界流束生成量γを算出する。本件の場合は、図１４に示すように、境界における流束（Ｊｆｌｕｘ）は、化学反応等により変化しないと仮定しているため、境界流束生成量γ＝０である。
【０１６０】
次に、第一領域側境界流束算出部７２は、上記式（４６）を用いて第一領域の境界近傍の流束ｆ−（λＣＬｋ＋１）を算出する。なお、本件では、α＝１であり、β＝１である。
【０１６１】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２より算出された流束ｆ−（λＣＬｋ＋１）を用いて上記式（５２）を解き、電界質内の含水量λＣＬｋ＋１を求める。
【０１６２】
そして、第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた電界質内の含水量λＣＬｋ＋１の値を用いて、第一領域の境界近傍における電界質内の含水量λの値λ−ｋ＋１を算出する。
【０１６３】
以上により、第一領域内の所定時間Δｔ経過後における電界質内の含水量λの値を全て求めることができる。
【０１６４】
次に、第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（４９）により境界における電界質内の含水量λの境界変化量ｃを算出する。本件では、境界におけるλの値は化学反応等により変化しないと仮定しているため、境界変化量ｃ＝０である。
【０１６５】
第二領域側境界物質量算出部８２は、上記式（５３）を用いて第二領域の境界近傍の電界質内の含水量λＣＬｋ＋１を算出する。なお、本件では、ａ＝１でありｂ＝１である。
【０１６６】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２より算出された電界質内の含水量λＣＬｋ＋１を用いて上記式（５１）を解き、電界質内の含水量λＧＤＬｋ＋１を求める。
【０１６７】
そして、第二領域側境界流束算出部８４は、第二領域内物理量解析部８３により求められた電界質内の含水量λＧＤＬｋ＋１の値を用いて、第二領域の境界近傍における電界質内の含水量λの流束ｇ＋（λＧＤＬｋ＋１）を算出する。
【０１６８】
以上により、第二領域内の所定時間Δｔ経過後における電界質内の含水量λｗの値を全て求めることができる。
【実施例２】
【０１６９】
『本発明における薄膜近似を用いた高分子形燃料電池の解析』
本実施例２では、上述した第二実施形態に基づく、本発明における薄膜近似を用いた高分子形燃料電池の解析手法について説明する。
【０１７０】
本実施例２における解析場は、図１５に示すように、電解膜、触媒層および多孔質層の３つの領域からなり、中間に位置する触媒層における物理量を薄膜近似して処理するものである。
【０１７１】
（１）各層における輸送方程式について
各層における輸送方程式について、図１５を参照しながら説明する。図１５に示すように、電解膜および多孔質層における液水ｓ、水蒸気質量分率Ｙｗおよび電解質内の含水量λの輸送方程式は、上述した実施例１におけるものと同一である。
【０１７２】
具体的には、電解膜内では、液水ｓ＝０、水蒸気質量分率Ｙｗ＝０であるとともに、電解質内の含水量λの輸送方程式は上記式（２５）により表される。また、多孔質層内では、液水ｓの輸送方程式は上記式（１８）により表され、水蒸気質量分率Ｙｗの輸送方程式は上記式（２３）により表されるとともに、電解質内の含水量λ＝０である。
【０１７３】
なお、以下の説明においては、上記式（２５）は、基礎方程式形で表された、以下の式（５４）として表した。
・・・（５４）
同様に、上記式（１８）は、以下の式（５５）とした。
・・・（５５）
さらに、上記式（２３）は、以下の式（５６）とした。
・・・（５６）
【０１７４】
また、本実施例２では、触媒層内の物理量の輸送は一定であると仮定して、その平均値ｓ、Ｙｗ、λおよび平均流束ｈ（ｓ）、ｈ（Ｙｗ）、ｈ（λ）として定義した。
【０１７５】
（２）各輸送方程式の解析について
本実施例２では、液水ｓ、水蒸気質量分率Ｙｗおよび電解質内の含水量λを同時に処理する。また、多孔質層を第一領域、触媒層を第二領域および電解膜を第三領域とした。
【０１７６】
電解膜から多孔質層への物理量の移動に対しては、ノイマン条件を与えた。ただし、電解膜における電解質内の含水量λは、実際には多孔質層内の水蒸気質量分率Ｙｗおよび液水ｓとして移動するが、すべてが多孔質層内の水蒸気質量分率Ｙｗとして移動するものとして仮定した。なお、この仮定下においては、多孔質層内の水蒸気質量分率Ｙｗが飽和量を超えた際には直ちに液水ｓに変換される関係にあり、電解質内の含水量λは間接的に液水ｓになるものであるため、物理的な関連性は保たれているものと考えられる。
【０１７７】
以上の仮定を踏まえると、第三領域側境界におけるノイマン条件は、以下の式（５７）
で表される。・・・（５７）
また、第一領域側境界におけるノイマン条件は、以下の式（５８）で表される。
・・・（５８）
【０１７８】
また、多孔質層から電解膜への物理量の移動に対しては、ディレクレ条件を与えた。ここで、触媒層内における電解質内の含水量λは、平衡状態に達していると考え、以下の式（５９）に示すように、液水ｓとａｗ（water activity）の関数で表されるものとした。
・・・（５９）
【０１７９】
第三領域側境界におけるディレクレ条件は、以下の式（６０）で表される。
・・・（６０）
また、第一領域境界においては、液水ｓおよび水蒸気質量分率Ｙｗのディレクレ条件が与えられ、液水ｓに関しては、以下の式（６１）で表される。
・・・（６１）
また、水蒸気質量分率Ｙｗに関しては、以下の式（６２）で表される。
・・・（６２）
【０１８０】
次に、上述した各式の離散化を行う。第一領域における電解質内の含水量λの輸送方程式である上記式（５４）の離散化をすると、以下の式（６３）として表される。
・・・（６３）
また、第三領域における液水ｓの輸送方程式である上記式（５５）の離散化をすると、以下の式（６４）として表される。
・・・（６４）
さらに、第三領域における水蒸気質量分率Ｙｗの輸送方程式である上記式（５６）を離散化すると、以下の式（６５）として表される。
・・・（６５）
【０１８１】
次に、各境界条件についても離散化を行う。第三領域側境界におけるノイマン条件の上記式（５７）を離散化すると、以下の式（６６）として表される。
・・・（６６）
また、本実施例２では、電解質内の含水量λは、すべてが多孔質層内の水蒸気質量分率Ｙｗとして移動するものと仮定したため、電解質内の含水量λと水蒸気質量分率Ｙｗとの間には、以下の式（６７）の関係が成り立つ。
・・・（６７）
また、第一領域側境界におけるノイマン条件の上記式（５８）の離散化を行うと、以下の式（６８）として表される。
・・・（６８）
【０１８２】
次に、第三領域側境界のディレクレ条件の上記式（６０）を離散化すると、以下の式（６９）として表される。
・・・（６９）
また、第三領域側境界における液水ｓのディレクレ条件の上記式（６１）の離散化をすると、以下の式（７０）として表される。
・・・（７０）
さらに、第一領域側境界における水蒸気質量分率Ｙｗのディレクレ条件の上記式（６２）を離散化すると、以下の式（７１）として表される。
・・・（７１）
【０１８３】
また、電解質内の含水量λを表した上記式（５９）を離散化すると、以下の式（７２）として表される。
・・・（７２）
【０１８４】
次に、以上の各式を用いて、第一領域および第三領域における各物理量の輸送方程式の解析について、特に、第二実施形態における第一領域内解析手段９および第三領域内解析手段１０の解析手順に従って、以下に説明する。
【０１８５】
第一領域内解析手段９では、まず、第一領域側境界流束生成量算出部９１が、境界流束生成量γ＋を算出する。本件では、境界流束生成量γ＋は、以下の式（７３）として算出される。
・・・（７３）
なお、陰極側の計算には、以下の式（７４）で表されるバトラー・ボルマー方程式を用いた。
・・・（７４）
ここで、電子２ｍｏｌ流れたときにＨ２Ｏは１ｍｏｌ発生するため、以下の式（７５）で表される水が発生する。
・・・（７５）
この式（７５）に、触媒層の厚みＬを掛けると、上記式（７３）の陰極側の値が導き出される。
【０１８６】
次に、第三領域側境界境界流束生成量算出部９２が、境界流束生成量γ−を算出する。本件では、境界流束生成量γ−＝０である。
【０１８７】
次に、第一領域側境界流束算出部９３は、上記式（６６）、上記式（６７）および上記式（６８）を用いて第一領域側境界の流束ｇ（Ｙｗ、ＧＤＬｋ＋１）を算出する。なお、本件におけるαおよびβを以下に示す。
【０１８８】
第一領域内物理量解析部９４は、第一領域側境界流束算出部７３より算出された流束ｇ（Ｙｗ、ＧＤＬｋ＋１）を用いて上記式（６５）を解き、多孔質層内の水蒸気質量分率Ｙｊｋ＋１を求める。
【０１８９】
そして、第一領域側境界物理量算出部９５は、第一領域内物理量解析部９４により求められた多孔質層内の水蒸気質量分率Ｙｊｋ＋１の値を用いて、第一領域側境界における水蒸気質量分率Ｙｗの値Ｙ＋ｋ＋１を算出する。
【０１９０】
以上により、第一領域内の所定時間Δｔ経過後における水蒸気質量分率Ｙｗの値を全て求めることができる。
【０１９１】
次に、第三領域内解析手段１０では、第三領域側境界変化量算出部１０１が、上記式（６０）により境界変化量ｃ−を算出する。本件では、境界変化量ｃ−＝０である。
【０１９２】
また、第一領域側境界変化量算出部１０２が、上記式（６１）および上記式（６２）により、境界変化量ｃｓ，＋および境界変化量ｃＹ，＋を算出する。本件では、境界変化量ｃｓ，＋＝０であり、境界変化量ｃＹ，＋＝０である。
【０１９３】
次に、第三領域側境界物質量算出部１０３は、上記式（６９）、上記式（７１）および上記式（７２）を用いて第三領域側境界の電界質内の含水量λ−ｋ＋１を算出する。なお、本件では、ａおよびｂに関する値はすべて１である。
【０１９４】
第三領域内物理量解析部１０４は、第三領域側境界物質量算出部１０３より算出された電界質内の含水量λ−ｋ＋１を用いて上記式（６３）を解き、電界質内の含水量λｉｋ＋１を求める。
【０１９５】
そして、第三領域側境界流束算出部１０５は、第三領域内物理量解析部１０４により求められた電界質内の含水量λｉｋ＋１の値を用いて、第三領域側境界における電界質内の含水量λの流束ｆ＋（λＧＤＬｋ＋１）を算出する。
【０１９６】
以上により、第三領域内の所定時間Δｔ経過後における電界質内の含水量λｗの値を全て求めることができる。
【０１９７】
なお、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法による解析対象は、前述した高分子燃料電池の問題に限定されるものではなく、他にも適用することができる。例えば、二次電池におけるイオン伝導問題や多孔質内における燃焼問題等に使用することができる。
【符号の説明】
【０１９８】
１ コンピュータ
１ａ，１ｂ 輸送方程式解析プログラム
２ 入力手段
３ 表示手段
４ 記憶手段
５ 演算手段
６ 物理量取得手段
７ 第一領域内解析手段
８ 第二領域内解析手段
９ 第一領域内解析手段
１０ 第三領域内解析手段
４１ プログラム記憶部
４２ 物理量記憶部
７１ 境界流束生成量算出部
７２ 第一領域側境界流束算出部
７３ 第一領域内物理量解析部
７４ 第一領域側境界物理量算出部
８１ 境界変化量算出部
８２ 第二領域側境界物理量算出部
８３ 第二領域内物理量解析部
８４ 第二領域側境界流束算出部
９１ 第一領域側境界流束生成量算出部
９２ 第三領域側境界流束生成量算出部
９３ 第一領域側境界流束算出部
９４ 第一領域内物理量解析部
９５ 第一領域側境界物理量算出部
１０１ 第三領域側境界変化量算出部
１０２ 第一領域側境界変化量算出部
１０３ 第三領域側境界物理量算出部
１０４ 第三領域内物理量解析部
１０５ 第三領域側境界流束算出部
【技術分野】
【０００１】
本発明は、それぞれ性質の異なる複数の領域を移動する所定の物理量の輸送方程式を解析する輸送方程式解析プログラムおよび輸送方程式解析方法に関するものである。
【背景技術】
【０００２】
近年、燃料電池は、電力発生に伴うＣＯ２排出のないクリーンな次世代エネルギーとして注目されている。この燃料電池は、多孔質層や電解膜等から構成されており、水素等の燃料剤と酸素等の酸化剤を化学反応させることにより電力を取り出すシステムである。従来、この燃料電池の発電効率の向上を目的として、燃料電池内の化学反応解析や流動解析等を行うプログラムの開発がなされている。
【０００３】
例えば、特開２００５−４４６０２号公報においては、燃料電池の物理的及び化学的特性の解析方法並びにプログラムが提案されている（特許文献１）。この特許文献１に記載のプログラムは、はじめに初期値の設定を行い、その後、酸素濃度、電流密度、水素濃度、温度分布、化学反応をそれぞれ解析し、これらのデータをもとに燃料電池全体の電位が釣り合うか否かの判別をし、さらに電流値が実験値と一致するか否かの判別を行うものであり、電気化学的拘束条件下の物質や熱の流れを容易に、かつ精度よく解析することができると記載されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００５−４４６０２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、特許文献１に記載された発明を含む、これまでの燃料電池の数値解析プログラムは以下のような問題点を有している。すなわち、
１．従来の数値解析プログラムは、発電セルの平面内における所定の物理量の分布を解析するものである。よって、それぞれ性質の異なる多層領域間の物理量の移動を正確に求めているわけではなく、多層領域間を移動する物理量の計算や当該物理量の計算に基づいた発電量の計算等の解析結果がその精度において不十分であった。
２．燃料電池などの多層領域間における物理量は、各層の境界における化学反応等により各種の物理量の変化を伴い、解析結果は前記境界における境界条件に強い影響を受ける。そのため、与えた境界条件や初期値によっては収束性が悪くなったり、解が得られない場合があった。
３．燃料電池の解析において、多孔質材と電解質膜との間にある触媒層は、他の層に比べて非常に薄く、解析精度は触媒層近傍における格子解像度に依存していた。そのため、格子を十分に小さくする必要があったが、その分、計算負荷が大きくなっていた。
【０００６】
本発明は、このような問題点を解決するためになされたものであって、解析精度の向上、解析速度の向上、解析負荷の軽減ならびに収束性の安定を図ることのできる輸送方程式解析プログラムおよび輸送方程式解析方法を提供することを目的としている。
【課題を解決するための手段】
【０００７】
本発明に係る輸送方程式解析方法を用いた輸送方程式解析プログラムは、それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析プログラムであって、前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出部と、前記境界流束生成量算出部により算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第二領域内解析手段は、前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出部と、前記境界変化量算出部により算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出部と、前記第二領域側境界物理量算出部により算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析部と、前記第二領域内物理量解析部により求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出部として機能させる。
【０００８】
すなわち、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法は、第一領域と第二領域との境界において、ディレクレ条件に基づく境界条件とノイマン条件に基づく境界条件の２つの境界条件式を用いるとともに、その２つの境界条件式を所定の順序で解くことによって、汎用的な輸送方程式解析プログラムおよび輸送方程式解析方法を構築したものである。
【０００９】
また、本発明に係る輸送方程式解析方法を用いた輸送方程式解析プログラムは、それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析プログラムであって、前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出部と、前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出部と、前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第三領域内解析手段は、前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出部と、前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出部と、前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出部と、前記第三領域側境界物理量算出部により算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析部と、前記第三領域内物理量解析部により求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出部として機能させる。
【００１０】
すなわち、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法は、第二領域が薄膜近似と捉えられる場合であって、前記第一領域側境界および前記第三領域側境界に、ディレクレ条件に基づく境界条件とノイマン条件に基づく境界条件の２つの境界条件式を用いるとともに、それらの合計４つの境界条件式を所定の順序で解くことによって、薄膜近似を用いた汎用的な輸送方程式解析プログラムおよび輸送方程式解析方法を構築したものである。
【発明の効果】
【００１１】
本発明によれば、解析精度の向上、解析速度の向上、解析負荷の軽減ならびに収束性の安定を図ることができる。
【図面の簡単な説明】
【００１２】
【図１】本発明に係る輸送方程式解析プログラムを備えたコンピュータの第一実施形態を示すブロック図である。
【図２】本第一実施形態における輸送方程式解析プログラムにより解析を行う第一領域および第二領域を示す概略図である。
【図３】本第一実施形態の輸送方程式解析プログラムおよび輸送方程式解析方法における処理の手順を示すフローチャート図である。
【図４】本第一実施形態における第一領域内解析手段の処理の手順を示すフローチャート図である。
【図５】本第一実施形態における第二領域内解析手段の処理の手順を示すフローチャート図である。
【図６】本第一実施形態における各構成の機能を示す機能ブロック図である。
【図７】本発明に係る輸送方程式解析プログラムを備えたコンピュータの第二実施形態を示すブロック図である。
【図８】本第二実施形態における輸送方程式解析プログラムにより解析を行う第一領域ないし第三領域を示す概略図である。
【図９】本第二実施形態の輸送方程式解析プログラムおよび輸送方程式解析方法における処理の手順を示すフローチャート図である。
【図１０】本第二実施形態における第一領域内解析手段の処理の手順を示すフローチャート図である。
【図１１】本第二実施形態における第三領域内解析手段の処理の手順を示すフローチャート図である。
【図１２】本第二実施形態における各構成の機能を示す機能ブロック図である。
【図１３】高分子形燃料電池の概略構成を示す概略構成図でる。
【図１４】本実施例１における電解膜、触媒層および多孔質層での液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λの輸送方程式および境界条件を示す解析場概略図である。
【図１５】本実施例２における電解膜、触媒層および多孔質層での液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λの輸送方程式および境界条件を示す解析場概略図である。
【発明を実施するための形態】
【００１３】
以下、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法の第一実施形態について図面を用いて説明する。図１は、本第一実施形態における輸送方程式解析プログラム１ａを備えたコンピュータの構成を示すブロック図である。
【００１４】
また、図２は、本第一実施形態で解析を行う第一領域および第二領域を示す概略図である。前記第一領域および前記第二領域は、境界を介して連結されており、それぞれ性質の異なる構造および／または物質からなる。また、前記各領域の間を移動する所定の物理量は、異なる輸送方程式により表される。例えば、第一領域における輸送方程式は、以下の式（１）により表される。
・・・式（１）
また、第二領域における輸送方程式は、以下の式（２）により表される。
・・・式（２）
ここで、上記式（１）における「φ」は第一領域における物理量であり、上記式（２）における「ψ」は第二領域における物理量である。なお、本発明において、上記「φ」と上記「ψ」は同じ物理量を表すものであり、上記式（１）および上記式（２）における物理量を同じ符号を用いて記載することも可能であるが、以下に記載される各式における物理的意味を明確にするため、第一領域における物理量と第二領域における物理量を別表記とした。
【００１５】
また、本発明における物理量は、例えば、水、水蒸気および熱等を示すものであり、各領域において輸送方程式により表されるもの全ての物理量である。
【００１６】
本第一実施形態におけるコンピュータ１Ａは、図１に示すように、主として、入力手段２、表示手段３、記憶手段４および演算処理手段５から構成されている。以下、各構成について詳細に説明する。
【００１７】
入力手段２は、テキストや数値、記号等を入力する操作キー、操作マウス等からなる。本第一実施形態においては、輸送方程式解析プログラム１ａを実行する際に用いられる初期条件や設定条件等の数値の入力操作等に用いることができるようになっている。
【００１８】
表示手段３は、画像やテキストデータを表示する液晶ディスプレーやＣＲＴディスプレー、タッチパネル等からなり、輸送方程式解析プログラム１ａにおけるユーザーインターフェースとして、入力手段２により入力された内容や解析結果等を表示できるようになっている。
【００１９】
記憶手段４は、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、ハードディスク、フラッシュメモリ等によって構成されており、各種のデータを記憶するとともに、演算処理手段５が演算を行う際のワーキングエリアとして機能するものである。
【００２０】
本第一実施形態において、記憶手段４は、図１に示すように、主として、輸送方程式解析プログラム１ａを記憶するプログラム記憶部４１と、空間的に離散化された各格子点における物理量を記憶させておく物理量記憶部４２とを有する。なお、本第一実施形態における物理量記憶部４２には、各格子点における物理量とともに、各格子点における流束が記憶されている。
【００２１】
演算処理手段５は、ＣＰＵ（Central Processing Unit）等から構成されており、記憶手段４にインストールされた第一実施形態の輸送方程式解析プログラム１ａを実行させることにより、図１に示すように、物理量取得手段６、第一領域内解析手段７および第二領域内解析手段８としてコンピュータ１Ａを機能させるようになっている。以下、各構成手段についてより詳細に説明する。
【００２２】
物理量取得手段６は、記憶手段４の物理量記憶部４２に記憶された所定の時間ｔの各格子点における物理量を取得するものである。また、本第一実施形態における物理量取得手段６は、各格子点における物理量とともに、各格子点の流束も取得している。なお、流束は取得した物理量から計算して求めることも可能である。
【００２３】
第一領域内解析手段７は、物理量取得手段６により取得された物理量の時間ｔにより、所定時間Δｔ後における第一領域内の物理量の数値解を求めるものである。例えば、物理量取得手段６により取得された物理量が初期値である場合、物理量取得手段６により取得された物理量の時間ｔは、ｔ＝０であり、第一領域内解析手段７は、その時間ｔ＝０から所定時間Δｔ後の物理量を求めるものである。
【００２４】
本第一実施形態における第一領域内解析手段７は、図１に示すように、境界流束生成量算出部７１と、第一領域側境界流束算出部７２と、第一領域内物理量解析部７３と、第一領域側境界物理量算出部７４とを有している。
【００２５】
境界流束生成量算出部７１は、以下の式（３）の流束境界条件式により定義された境界流束生成量γに係る関数γ（φ，ψ）を用いて前記境界流束生成量γの値を算出するものである。
・・・（３）
【００２６】
よって、境界流束生成量算出部７１では、前記境界流束生成量γの関数γ（φ，ψ）に、前記物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界流束生成量γを算出している。
【００２７】
なお、この境界流束生成量γは、境界における化学反応等により物理量が変化したことにより生じる流束の差を示す値であり、算出された境界流束生成量γの符号が正の場合には単位時間当たりの生成量を、符号が負の場合には単位時間当たりの消滅量を表す。
【００２８】
また、上記式（３）で示される流束境界条件式は、第一領域と第二領域との間のノイマン条件による境界条件を表す式であり、前記第一領域側から前記境界へ流入出する前記物理量の流束ｆ（φ）と前記第二領域側から前記境界へ流入出する前記物理量の流束ｇ（ψ）との連続性を表したものである。
【００２９】
ここで、「α」は第一領域における係数であり、「β」は第二領域における係数である。例えば、物理量が液水であって第一領域が多孔質層である場合、液水は多孔質層内の孔の中しか通過することができない。よって、第一領域における流束は、流束ｎｆ（φ）と第一領域の多孔率（孔の体積／全体の体積）εとの積で表される。よって、この場合はα＝εとなる。
【００３０】
第一領域側境界流束算出部７２は、前記境界流束生成量算出部７１により算出された境界流束生成量γ、および前記物理量取得手段により取得された第二領域における境界近傍の流束ｇ＋（ψｋ）を、流束境界条件式に代入して所定時間Δｔ経過後の第一領域における境界近傍の流束ｆ−（φｋ＋１）を算出するものである。よって、第一領域側境界流束算出部７２における流束境界条件式は、所定時間Δｔ経過後の流束ｆ−（φｋ＋１）を算出するため、以下の式（４）に示すように時間方向に離散化されている。
・・・式（４）
【００３１】
ここで、「ｋ」は時間を表すものであって物理量取得手段６により取得された物理量の有する時間ｔに対応するものである。また、「ｋ＋１」は、時間ｔから所定時間Δｔ経過後の時間を表すものである。
【００３２】
また、「＋」は第二領域における境界近傍の値であることを示しており、具体的には、境界から第二領域側にある一つ目の格子点における値を示している。また、「−」は第一領域における境界近傍の値であることを示しており、具体的には、境界から第一領域側にある一つ目の格子点における値を示している。
【００３３】
よって、式（４）は、第二領域から境界へと流入する時間ｋにおける流束ｇ＋（ψｋ）が、境界において境界流束生成量γ分だけ増減され、所定時間Δｔ経過後の時間ｋ＋１には、第一領域における境界近傍へとｆ−（φｋ＋１）として流出していることを示している。
【００３４】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２により算出された境界近傍の流束αｎｆ−（φｋ＋１）を第一領域の輸送方程式に代入して所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求めるものである。よって、第一領域内物理量解析部７３における第一領域の輸送方程式は、所定時間Δｔ経過後の第一領域内における物理量の数値解を求めるために空間的に離散化すると、上記式（１）を以下の式（５）として表される。
・・・（５）
ここで、ｉは、第一領域内の格子点の位置を表すものであり、本第一実施形態においては、境界から２番目の格子点から、ｉ＝１，２，３，・・・と番号が割り当てられる。
【００３５】
本第一実施形態における第一領域内物理量解析部７３は、第一領域側境界流束算出部７２により算出した第一領域における境界近傍の流束ｆ−（φｋ＋１）を上記式（５）に代入して導かれる連立方程式を解くことで、第一領域内の各格子点における所定時間Δｔ経過後の物理量φｉｋ＋１を算出する。なお、上記式（５）による連立方程式を解く方法は特に限定されるものではなく、陰解法や陽解法等から適宜選択されるものである。
【００３６】
第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の物理量ａφ−ｋ＋１を算出するものである。
【００３７】
次に、第二領域内解析手段８について説明する。第二領域内解析手段８は、所定時間Δｔ経過後における第二領域内の物理量の数値解を求めるものである。本第一実施形態における第二領域内解析手段８は、図１に示すように、境界変化量算出部８１と、第二領域側境界物理量算出部８２と、第二領域内物理量解析部８３と、第二領域側境界流束算出部８４とを有している。
【００３８】
境界変化量算出部８１は、以下の式（６）の物理量境界条件式により定義された境界変化量ｃに係る関数ｃ（φ，ψ）を用いて前記境界変化量ｃの値を算出するものである。
・・・（６）
【００３９】
よって、境界変化量算出部８１では、前記境界変化量ｃの関数ｃ（φ，ψ）に、前記物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界変化量ｃを算出している。なお、この境界変化量ｃは、境界における流束の伝達率により生じる第一領域側の物理量と第二領域側の物理量との差を示す値である。
【００４０】
また、上記式（６）で示される物理量境界条件式は、第一領域と第二領域との間のディレクレ条件による境界条件を表す式であり、境界における第一領域側の物理量φと第二領域側の物理量ψとの連続性を示している。
【００４１】
第二領域側境界物理量算出部８２と境界変化量算出部８１とにより算出された境界変化量ｃおよび第一領域における境界近傍の物理量φ−ｋを物理量境界条件式に代入して、所定時間Δｔ経過後の第二領域における境界近傍の物理量ψ＋ｋ＋１を算出するものである。よって、第二領域側境界物理量算出部８２における物理量境界条件式は、所定時間Δｔ経過後の物理量ｂ（ψｋ）ψ＋ｋ＋１を算出するため、以下の式（７）に示すように時間方向に離散化されている。
・・・（７）
【００４２】
ここで、「ａ」は第一領域における係数であり第一領域の物理量φの関数として与えられている。また、「ｂ」は第二領域における係数であり第二領域の物理量ψの関数として与えられる。
【００４３】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２により算出された境界近傍の物理量ｂ（ψｋ）ψ＋ｋ＋１を第二領域の輸送方程式に代入して所定時間Δｔ経過後の第二領域内における物理量ψｊｋ＋１の数値解を求めるものである。よって、第二領域内物理量解析部８３における第二領域の輸送方程式は、所定時間Δｔ経過後の第二領域内における物理量ψｊｋ＋１の数値解を求めるために空間的に離散化すると、上記式（２）を以下の式（８）として表される。
・・・（８）
ここで、ｊは、第二領域内の格子点の位置を表すものであり、本第一実施形態においては、境界から２番目の格子点から、ｊ＝１，２，３，・・・と番号が割り当てられる。
【００４４】
本第一実施形態における第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２により算出した第二領域における境界近傍の物理量ψ＋ｋ＋１を上記式（８）に代入して導かれる連立方程式を解くことで、第二領域内の各格子点における所定時間Δｔ経過後の物理量ψｊｋ＋１を算出している。なお、上記式（８）による連立方程式を解く方法は特に限定されるものではなく、陰解法や陽解法等から適宜選択されるものである。
【００４５】
第二領域側境界流束算出部８４は、第二領域内物理量解析部８３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の流束ｇ＋（ψｋ＋１）を算出するものである。
【００４６】
次に、本第一実施形態の輸送方程式解析プログラム１ａにおける各構成の作用および輸送方程式解析方法について、図３ないし図５に示すフローチャート図および図６示す各構成の機能を示す機能ブロック図を用いて説明する。
【００４７】
図３および図６に示すように、本第一実施形態の輸送方程式解析プログラム１ａおよび輸送方程式解析方法では、まず、物理量取得手段６が、物理量記憶部４２に記憶された各格子点における物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を取得する（ステップＳ１）。続いて、サブルーチン化された第一領域内解析手段７が、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求める（ステップＳ２）。続いて、サブルーチン化された第二領域内解析手段８は、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ後の第二領域内における物理量ψｊｋ＋１の数値解を求める（ステップＳ３）。そして、時間ｔが所望する時間Ｔに達したか否かが判別される（ステップＳ４）。
【００４８】
なお、本第一実施形態では、図３に示すように、第一領域内解析手段７の処理後に第二領域内解析手段８の処理を行っているが、どちらを先に処理してもよく、並列に処理してもよい。また、本第一実施形態では、第一領域内解析手段７および第二領域内解析手段８をサブルーチン化して解析処理を実行しているが、サブルーチン化せずに物理量取得手段６と一連で処理するようにしてもよい。
【００４９】
以下、第一領域内解析手段７および第二領域内解析手段８について、より詳細に説明する。
【００５０】
図４に示すように、ステップＳ２における第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（３）により表される流束境界条件式において定義される境界流束生成量γに係る関数γ（φ，ψ）に、物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界流束生成量γの算出を行う（ステップＳ２１）。算出された境界流束生成量γは、図６に示すように、次のステップの第一領域側境界流束算出部７２に送られる。
【００５１】
第一領域側境界流束算出部７２は、時間方向に離散化された上記式（４）により表される流束境界条件式に、境界流束生成量算出部７１により算出された境界流束生成量γ、および物理量取得手段６により取得された第二領域における境界近傍の流束ｇ＋（ψｋ）を代入して、所定時間Δｔ経過後の第一領域における境界近傍の流束ｆ−（φｋ＋１）を算出する（ステップＳ２２）。算出された第一領域における境界近傍の流束ｆ−（φｋ＋１）は、図６に示すように、次のステップの第一領域内物理量解析部７３に送られるとともに、物理量記憶部４２に送られる。
【００５２】
物理量記憶部４２に送られた所定時間Δｔ経過後の第一領域における境界近傍の流束ｆ−（φｋ＋１）は、第一領域における境界近傍の格子点における流束値として置き換えられる。このように、物理量記憶部４２内の物理量を所定時間Δｔ経過後の物理量に置き換えることにより、時間発展計算を行うことができる。
【００５３】
次に、第一領域内物理量解析部７３は、時間方向に離散化された上記式（５）により表される第一領域の輸送方程式に、第一領域側境界流束算出部７２により算出された境界近傍の流束ｆ−（φｋ＋１）を代入して導き出される連立方程式を解くことで、第一領域内の各格子点における所定時間Δｔ経過後の物理量φｉｋ＋１の算出を行う（ステップＳ２３）。算出された第一領域における物理量φｉｋ＋１は、図６に示すように、次のステップの第一領域側境界物理量算出部７４に送られるとともに、物理量記憶部４２に送られる。
【００５４】
第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の物理量φ−ｋ＋１を算出する（ステップＳ２４）。そして、算出された第一領域における境界近傍の物理量φ−ｋ＋１は、図６に示すように、物理量記憶部４２に送られる。以上で、第一領域内解析手段７によるステップＳ２は終了し、次の第二領域内解析手段８によるステップＳ３へと進む。
【００５５】
図５に示すように、ステップＳ３における第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（６）の物理量境界条件式により定義された境界変化量ｃに係る関数ｃ（φ，ψ）に、物理量取得手段６により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界変化量ｃの算出を行う（ステップＳ３１）。算出された境界変化量ｃは、図６に示すように、次のステップの第二領域側物理量算出部８２に送られる。
【００５６】
第二領域側境界物理量算出部８２は、時間方向に離散化された上記式（７）により表される物理量境界条件式に、境界変化量算出部８１により算出された境界変化量ｃ、および物理量取得手段６により取得された第一領域における境界近傍の物理量φ−ｋを代入して、所定時間Δｔ経過後の第二領域における境界近傍の物理量ψ＋ｋ＋１を算出する（ステップＳ３２）。算出された第二領域における境界近傍の物理量ψ＋ｋ＋１は、図６に示すように、次のステップの第二領域内物理量解析部８３に送られるとともに、物理量記憶部４２に送られる。
【００５７】
次に、第二領域内物理量解析部８３は、時間方向に離散化された上記式（８）により表される第二領域の輸送方程式に、第二領域側境界物理量算出部８２により算出された境界近傍の物理量ψ＋ｋ＋１を代入して導き出された連立方程式を解くことによって第二領域内の各格子点における所定時間Δｔ経過後の物理量ψｊｋ＋１の算出を行う（ステップＳ３３）。算出された第二領域における物理量ψｊｋ＋１は、図６に示すように、次のステップの第二領域側境界流束算出部８４に送られるとともに、物理量記憶部４２に送られる。
【００５８】
第一領域側境界流束算出部８４は、第二領域内物理量解析部８３により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の流束ｇ＋（ψｋ＋１）を算出する（ステップＳ３４）。そして、算出された第二領域における境界近傍の流束ｇ＋（ψｋ＋１）は、図６に示すように、物理量記憶部４２に送られる。以上で、第二領域内解析手段８によるステップＳ３は終了し、次のステップＳ４へと進む。
【００５９】
ステップＳ４では、時間ｔが所望する時間Ｔに達したか否かを行う。時間ｔが所望する時間Ｔに達していない場合（ステップＳ４：ＮＯ）、時間ｔ＝ｔ＋ΔｔとしてステップＳ１へもどり、ステップ１からステップＳ３の処理を繰り返し行う。このように、繰り返し処理を行うことにより第一領域内および第二領域内の各格子点における物理量φ，ψの時間発展計算が行われる。そして、時間ｔが所望する時間Ｔに達した場合（ステップＳ４：ＹＥＳ）、処理を終了させる。
【００６０】
以上のような本第一実施形態の輸送方程式解析プログラム１ａおよび輸送方程式解析方法によれば、以下のような効果を得ることができる。
１．多層領域間を移動する物理量の数値解析を行うことができる。
２．流束境界条件式および物理量境界条件式の２つの境界条件を算出する式を用いることにより、適当な各領域間の境界における物理量を求めることができるため、各領域における輸送方程式の数値解の精度を向上させるとともに、収束性も向上する。
３．第一領域側と第二領域側との解析を並列処理することができるため、処理時間を向上させることができる。
４．２つの輸送方程式の連成問題であれば物理量の種類については特に限定がなく、応用範囲が広い。
【００６１】
次に、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法の第二実施形態について説明する。なお、本第二実施形態のうち、前述した第一実施形態の構成と同一若しくは相当する構成については同一の符号を付して再度の説明を省略する。
【００６２】
図７は、本第二実施形態における輸送方程式解析プログラム１ｂを備えたコンピュータの構成を示すブロック図である。
【００６３】
また、図８は、本第二実施形態で解析を行う第一領域、第二領域および第三領域を示す概略図である。図８に示すように、本第二実施形態における第一領域、第二領域および第三領域は、第一領域と第二領域とがこれらの間の第一領域側境界を介して連結されており、第二領域と第三領域とがこれらの間の第三領域側境界を介して連結されている。
【００６４】
また、第一領域および第三領域は、それぞれ性質の異なる構造および／または物質からなり、前記各領域の間を移動する所定の物理量は、異なる輸送方程式により表されるものである。具体的には、第一領域における輸送方程式は、上記式（１）により表されるととともに、第三領域における輸送方程式は、上記式（２）により表される。
【００６５】
さらに、本第二実施形態における第二領域は薄膜状に形成されている。そのため、第二領域内における物理量および流束は、空間的にほぼ一定であると見なすことができる。よって、本第二実施形態では、第二領域内の物理量およびその物理量の流束を、それぞれθおよびｈ（θ）と近似した。本発明では、以上のような第二領域における物理量等の近似を薄膜近似と定めている。
【００６６】
以下、本第二実施形態の各構成について、より詳細に説明する。
【００６７】
本第二実施形態におけるコンピュータ１Ｂは、図７に示すように、主として、入力手段２、表示手段３、記憶手段４および演算処理手段５から構成されている。
【００６８】
このうち本第二実施形態における演算処理手段５は、記憶手段４にインストールされた第二実施形態の輸送方程式解析プログラム１ｂを実行させることにより、図７に示すように、物理量取得手段６、第一領域内解析手段９および第三領域内解析手段１０としてコンピュータ１Ｂを機能させるようになっている。
【００６９】
本第二実施形態における第一領域内解析手段９は、上述した第一実施形態における第一領域内解析手段７と同様に、所定時間Δｔ経過後における第一領域内の物理量の数値解を求めるものである。
【００７０】
また、本第二実施形態における第一領域内解析手段９は、図７に示すように、第一領域側境界流束生成量算出部９１と、第三領域側境界流束生成量算出部９２と、第一領域側境界流束算出部９３と、第一領域内物理量解析部９４と、第一領域側境界物理量算出部９５とを有している。つまり、上述した第一実施形態では、境界が１つしかないため境界の境界流束生成量の算出は、境界流束生成量算出部７１による１つの処理で行うことができるのに対し、本第二実施形態では、境界が、第一領域側境界と第三領域側境界の２つあるため、各境界における流束の算出についても第一領域側境界流束生成量算出部９１と第三領域側境界流束生成量算出部９２とによる２つの処理を要する。
【００７１】
第一領域側境界流束生成量算出部９１は、以下の式（９）の流束境界条件式により定義された境界流束生成量γ−に係る関数γ−（φ，θ）を算出するものであり、前記関数γ−（φ，θ）に取得された物理量φおよび物理量θを代入して第一領域側境界における物理量の境界流束生成量γ−を算出している。
・・・（９）
【００７２】
また、上記式（９）で示される流束境界条件式は、第一領域と第二領域との間のノイマン条件による境界条件を表す式であり、第一領域側から第一領域側境界へ流入出する物理量の流束ｆ（φ）と第二領域側から第一領域側境界へ流入出する前記物理量の流束ｈ（θ）との連続性を示している。ここで、「−」は、第一領域側境界の値であることを示している。
【００７３】
次に、第三領域側境界流束生成量算出部９２は、以下の式（１０）の流束境界条件式により定義された境界流束生成量γ＋に係る関数γ＋（θ，ψ）を算出するものであり、前記関数γ＋（θ，ψ）に物理量取得手段６により取得された物理量ψおよびθを代入して第三領域側境界における物理量の境界流束生成量γ＋を算出している。
・・・（１０）
【００７４】
また、上記式（１０）で示される流束境界条件式は、第三領域と第二領域との間のノイマン条件による境界条件を表す式であり、第三領域側から第三領域側境界へ流入出する物理量の流束ｇ（ψ）と第二領域側から第三領域側境界へ流入出する前記物理量の流束ｈ（θ）との連続性を示している。ここで、「＋」は、第三領域側境界の値であることを示している。
【００７５】
第一領域側境界流束算出部９３は、第一領域側境界流束生成量算出部９１により算出された境界流束生成量γ−、第三領域側境界流束生成量算出部９２により算出された境界流束生成量γ＋、および前記物理量取得手段６により取得された第三領域側境界の流束ｇ＋（ψｋ）を、前記流束境界条件式に代入して所定時間Δｔ経過後の第一領域側境界の流束ｆ−（φｋ＋１）を算出するものである。そのため、第一領域側境界流束算出部９３では、時間方向に離散化された、以下の式（１１）および式（１２）からなる流束境界条件式が用いられる。
・・・（１１）
・・・（１２）
【００７６】
具体的には、まず、上記式（１１）に第三領域側境界流束生成量算出部９２により算出された境界流束生成量γ＋と、第三領域境界の流束ｇ＋（ψｋ）を代入して、所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を算出する。そして、上記式（１２）に、第一領域側境界流束生成量算出部９１により算出された境界流束生成量γ−と式（１１）により算出された所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を代入して、所定時間Δｔ経過後の第一領域側境界の流束ｆ−（φｋ＋１）を算出する。
【００７７】
本第二実施形態における第一領域内物理量解析部９４は、第一実施形態における第一領域内物理量解析部７３と同様の機能を有しており、第一領域側境界流束算出部９３により算出された第一領域側境界の流束ｆ−（φｋ＋１）を、上記式（５）で表される第一領域の輸送方程式に代入して所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求めるものである。
【００７８】
また、本第二実施形態における第一領域側境界物理量算出部９５は、第一実施形態における第一領域側境界物理量算出部７４と同様の機能を有しており、第一領域内物理量解析部９４により求められた数値解から所定時間Δｔ経過後の第一領域側境界の物理量φ−ｋ＋１を算出するものである。
【００７９】
次に、第三領域内解析手段１０について説明する。本第二実施形態における第三領域内解析手段１０は、第一実施形態における第二領域内解析手段８と同様に、所定時間Δｔ経過後における第三領域内の物理量の数値解を求めるものである。本第二実施形態における第二領域内解析手段１０は、図７に示すように、第三領域側境界変化量算出部１０１と、第一領域側境界変化量算出部１０２と、第三領域側境界物理量算出部１０３と、第三領域内物理量解析部１０４と、第三領域側境界流束算出部１０５とを有している。
【００８０】
第三領域側境界変化量算出部１０１は、以下の式（１３）の流束境界条件式により定義された境界変化量ｃ＋に係る関数ｃ＋（θ，ψ）を算出するものであり、前記関数ｃ＋（θ，ψ）に前記物理量取得手段６により取得された物理量ψおよび物理量θを代入して第三領域側境界における物理量の境界変化量ｃ＋を算出している。
・・・（１３）
【００８１】
また、上記式（１３）で示される流束境界条件式は、第三領域と第二領域との間のディレクレ条件による境界条件を表す式であり、第三領域側境界における第三領域側の物理量ψと第二領域側の物理量θとの連続性を示している。
【００８２】
第一領域側境界変化量算出部１０２は、以下の式（１４）の流束境界条件式により定義された境界変化量ｃ−に係る関数ｃ−（φ，θ）を算出するものであり、前記関数ｃ−（φ，θ）に前記物理量取得手段６により取得された物理量φおよび物理量θを代入して第一領域側境界における物理量の境界変化量ｃ−を算出している。
・・・（１４）
【００８３】
また、上記式（１４）で示される流束境界条件式は、第一領域と第二領域との間のディレクレ条件による境界条件を表す式であり、第一領域側境界における第一領域側の物理量φと第二領域側の物理量θとの連続性を示している。
【００８４】
第三領域側境界物理量出部１０３は、第三領域側境界変化量算出部１０１により算出された境界変化量ｃ＋、第一領域側境界変化量算出部１０２により算出された境界変化量ｃ−、および前記物理量取得手段６により取得された第一領域境界の物理量φ−を、物理量境界条件式に代入して所定時間Δｔ経過後の第三領域側境界の物理量ψ＋ｋ＋１を算出するものである。そのため、第三領域側境界物理量算出部１０３では、時間方向に離散化された、以下の式（１５）および式（１６）からなる物理量境界条件式が用いられる。
・・・（１５）
・・・（１６）
【００８５】
具体的には、まず、上記式（１５）に第一領域側境界変化量算出部１０２により算出された境界変化量ｃ−と、第一領域境界の物理量φ−を代入して、所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を算出する。そして、上記式（１６）に、第三領域側境界変化量算出部１０１により算出された境界変化量ｃ＋と、式（１５）により算出された所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を代入して、所定時間Δｔ経過後の第三領域側境界の物理量ψ−ｋ＋１を算出する。
【００８６】
本第二実施形態における第三領域内物理量解析部１０４は、第一実施形態における第二領域内物理量解析部８３と同様の機能を有しており、第三領域側境界物理量算出部１０３により算出された第三領域側境界の物理量ψ−ｋ＋１を、上記式（８）で表される第三領域の輸送方程式に代入して所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求めるものである。
【００８７】
また、本第二実施形態における第三領域側境界流束算出部１０５は、第一実施形態における第二領域側境界流束算出部８４と同様の機能を有しており、第三領域内物理量解析部１０４により求められた数値解から所定時間Δｔ経過後の第三領域側境界の物理量の流束ｇ＋（ψｋ＋１）を算出するものである。
【００８８】
次に、本第二実施形態の輸送方程式解析プログラム１ｂにおける各構成の作用および輸送方程式解析方法について、図９ないし図１１に示すフローチャート図および図１２に示す各構成の機能を示す機能ブロック図を用いて説明する。
【００８９】
図９および図１２に示すように、本第二実施形態の輸送方程式解析プログラム１ｂおよび輸送方程式解析方法では、まず、物理量取得手段６が、物理量記憶部４２に記憶された各格子点における物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を取得する（ステップＳＳ１）。続いて、サブルーチン化された第一領域内解析手段９が、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ経過後の第一領域内における物理量φｉｋ＋１の数値解を求める（ステップＳＳ２）。続いて、サブルーチン化された第三領域内解析手段１０は、取得された物理量φ，ψおよび流束ｆ（φ），ｇ（ψ）を用いて所定時間Δｔ経過後の第三領域内における物理量ψｊｋ＋１の数値解を求める（ステップＳＳ３）。そして、時間ｔが所望する時間Ｔに達したか否かが判別される（ステップＳＳ４）。
【００９０】
なお、第一領域内解析手段９の処理と第三領域内解析手段１０の処理とはどちらを先に処理してもよく、並列に処理してもよい。また、本第二実施形態では、第一領域内解析手段９および第三領域内解析手段１０をサブルーチンとせずに物理量取得手段６と一連で処理するようにしてもよい。
【００９１】
図１０に示すように、ステップＳＳ２における第一領域内解析手段９では、まず、第一領域側境界流束生成量算出部９１が、上記式（９）の第一領域側境界における流束境界条件式により定義された境界流束生成量γ−に係る関数γ−（φ，θ）に、物理量取得手段６により取得された物理量φおよび物理量θを代入して、第一領域側境界における物理量の境界流束生成量γ−の算出を行う（ステップＳＳ２１）。算出された境界流束生成量γ−は、図１２に示すように、第一領域側境界流束算出部９３に送られる。
【００９２】
また、第三領域側境界流束生成量算出部９２が、上記式（１０）により表される第三領域側境界における流束境界条件式により定義され境界流束生成量γ＋に係る関数γ＋（θ，ψ）に、物理量取得手段６で取得された物理量ψおよびθを代入して、第一領域側境界における物理量の境界流束生成量γ＋の算出を行う（ステップＳＳ２２）。算出された境界流束生成量γ＋は、図１２に示すように、第一領域側境界流束算出部９３に送られる。
【００９３】
第一領域側境界流束算出部９３は、時間方向に離散化された上記式（１１）により表される第三領域側境界における流束境界条件式に、第三領域側境界流束生成量算出部９２により算出された境界流束生成量γ＋、および物理量取得手段６により取得された第三領域側境界の流束ｇ＋（ψｋ）を代入して、所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を算出するとともに、上記式（１２）に、第一領域側境界流束生成量算出部９１により算出された境界流束生成量γ−と式（１１）により算出された所定時間Δｔ経過後の第二領域内の流束ｈ（θｋ＋１）を代入して、所定時間Δｔ経過後の第一領域側境界の流束ｆ−（φｋ＋１）を算出する（ステップＳＳ２３）。算出された第一領域側境界の流束ｆ−（φｋ＋１）は、図１２に示すように、次のステップの第一領域内物理量解析部９４に送られるとともに、物理量記憶部４２に送られる。
【００９４】
次に、第一領域内物理量解析部９４は、時間方向に離散化された上記式（５）により表される第一領域の輸送方程式に、第一領域側境界流束算出部９３により算出された第一領域側境界の流束ｆ−（φｋ＋１）を代入して導き出される連立方程式を解くことにより、第一領域内の各格子点における所定時間Δｔ経過後の物理量φｉｋ＋１の算出を行う（ステップＳＳ２４）。算出された第一領域における物理量φｉｋ＋１は、図１２に示すように、次のステップの第一領域側境界物理量算出部９５に送られるとともに、物理量記憶部４２に送られる。
【００９５】
第一領域側境界物理量算出部９５は、第一領域内物理量解析部９４により求められた数値解から所定時間Δｔ経過後の第一領域における境界近傍の物理量φ−ｋ＋１を算出する（ステップＳＳ２５）。そして、算出された第一領域における境界近傍の物理量φ−ｋ＋１は、図１２に示すように、物理量記憶部４２に送られる。
【００９６】
図１１に示すように、ステップＳＳ３における第二領域内解析手段１０では、まず、第三領域側境界変化量算出部１０１が、上記式（１３）の第三領域側境界における物理量境界条件式により定義された境界変化量ｃ＋（θ，ψ）に係る関数に、物理量取得手段６により取得された物理量ψおよび物理量θを代入して、第三領域側境界における物理量の境界変化量ｃ＋の算出を行う（ステップＳＳ３１）。算出された境界変化量ｃ＋は、図１１に示すように、第三領域側物理量算出部１０３に送られる。
【００９７】
次に、第一領域側境界変化量算出部１０２が、上記式（１４）の第一領域側境界における物理量境界条件式により定義された境界変化量ｃ−（φ，θ）に係る関数に、物理量取得手段６により取得された物理量φおよび物理量θを代入して、第一領域側境界における物理量の境界境界変化量ｃ−の算出を行う（ステップＳＳ３２）。算出された境界変化量ｃ−は、図１１に示すように、第三領域側物理量算出部１０３に送られる。
【００９８】
第三領域側境界物理量算出部１０３は、時間方向に離散化された上記式（１５）により表される第一領域側境界における流束境界条件式に、第一領域側境界変化量算出部１０２により算出された境界変化量ｃ−、および物理量取得手段６により取得された第一領域側境界の物理量φｋを代入して、所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を算出するとともに、上記式（１６）に、第三領域側境界変化量算出部１０１により算出された境界変化量ｃ＋と、式（１５）により算出された所定時間Δｔ経過後の第二領域内の物理量θｋ＋１を代入して、所定時間Δｔ経過後の第三領域側境界の物理量ψ＋ｋ＋１を算出する（ステップＳＳ２３）。算出された第三領域側境界の物理量ψ＋ｋ＋１は、図１２に示すように、次のステップの第三領域内物理量解析部１０４に送られるとともに、物理量記憶部４２に送られる。
【００９９】
次に、第三領域内物理量解析部１０４は、時間方向に離散化された上記式（８）により表される第三領域の輸送方程式に、第三領域側境界物理量算出部１０３により算出された第三領域側境界における物理量ψ＋ｋ＋１を代入して導き出される連立方程式を解くことにより、第三領域内の各格子点における所定時間Δｔ経過後の物理量ψｊｋ＋１の算出を行う（ステップＳＳ３４）。算出された第三領域側境界における物理量ψｊｋ＋１は、図１２に示すように、次のステップの第三領域側境界流束算出部１０５に送られるとともに、物理量記憶部４２に送られる。
【０１００】
第三領域側境界流束算出部１０５は、第三領域内物理量解析部１０４により求められた数値解から所定時間Δｔ経過後の第三領域側境界における流束ｇ＋（ψｋ＋１）を算出する（ステップＳＳ３５）。そして、算出された第三領域側境界の流束ｇ＋（ψｋ＋１）は、図１２に示すように、物理量記憶部４２に送られる。
【０１０１】
ステップＳＳ４では、時間ｔが所望する時間Ｔに達したか否かを行う。時間ｔが所望する時間Ｔに達していない場合（ステッＳプＳ４：ＮＯ）、時間ｔ＝ｔ＋Δｔとし、ステップＳＳ１からステップＳＳ３の処理を繰り返し行う。そして、時間ｔが所望する時間Ｔに達した場合（ステップＳＳ４：ＹＥＳ）、処理を終了させる。
【０１０２】
以上のような本第二実施形態の輸送方程式解析プログラム１ｂおよび輸送方程式解析方法によれば、薄膜近似した第二領域にける格子解像度を高くする必要が無くなるため、計算負荷を小さくすることができるという効果が得られる。
【実施例１】
【０１０３】
『本発明を用いた高分子形燃料電池の解析』
本実施例１では、上述した第一実施形態の解析方法を用いて、高分子形燃料電池内における各種物理量の輸送方程式および各層の境界における境界条件式を例示しつつ、解析について説明する。
【０１０４】
高分子形燃料電池は、図１３に示すように、電解膜と、この電界膜を中心に陰極側および陽極側にそれぞれ触媒層と多孔質層とが多層状に連結されている。また、高分子形燃料電池では、水素や酸素などの作動ガス、生成物である水または水蒸気等が混合された状態で前記各層内を移動する、いわゆる複合流動場が構成されている。
【０１０５】
そこで、本実施例１では、各層における物理量を液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λに分けて、それぞれについて輸送方程式が設定されている。また、液水ｓおよび水蒸気質量分率Ｙｗは直接相変化が行われるとともに、水蒸気質量分率Ｙｗおよび電界質内の含水量λは直接相変化が行われるものとし、一方、液水ｓおよび電界質内の含水量λは直接相変化が行われないものとした。
【０１０６】
（１）各層における輸送方程式について
各層における輸送方程式について、図１４を参照しながら説明する。図１４に示すように、以下の式（１７）は、触媒層内における液水ｓの輸送方程式である。
・・・（１７）
また、以下の式（１８）は、多孔質層内における液水ｓの輸送方程式である。
・・・（１８）
【０１０７】
ここで、上記式（１７）および上記式（１８）内の各係数は、以下の式（１９）ないし式（２１）の関係を有している。
・・・（１９）
・・・（２０）
・・・（２１）
【０１０８】
また、εは多効率、ρｌは水の密度、ρｇはガスの密度、ｕ気体の速度、Ｄｃ拡散係数、κ透過率、νは動粘性計数、ｇは重力加速度、Ｍｌはモル質量、ｋｃは凝縮係数、ｋｅは蒸発係数、Ｘｗはモル分率、Ｒは気体定数、Ｔ温度、Ｐ圧力、Ｐｓａｔは飽和水蒸気量、σは表面張力、θｃは液水ｓの接触角およびμは粘性係数を示している。また、ｍ１は液水ｓと水蒸気質量分率Ｙｗの間の相変化量である。
【０１０９】
なお、電界膜は固体であり液水ｓは入り込めないため、図１４に示すように、電界膜内の液水ｓ＝０である。
【０１１０】
次に、水蒸気質量分率Ｙｗの輸送方程式について説明する。以下の式（２２）は、触媒層内における水蒸気質量分率Ｙｗの輸送方程式である。
・・・（２２）
また、以下の式（２３）は、多孔質層内における水蒸気質量分率Ｙｗの輸送方程式である。
・・・（２３）
【０１１１】
また、上記式（２２）および上記式（２３）内の各係数は、以下の式（２４）の関係を有している。
・・・（２４）
【０１１２】
ここで、Ｙｉは化学種ｉの質量分率、Ｄｉ，ｅｆｆは参照拡散係数、Ｄｉ０は拡散係数、Ｐ０は基準圧力、Ｔ０は基準温度、ωは化学反応による水または水蒸気の発生量、ｍ２は水蒸気質量分率Ｙｗと電界質内の含水量λと間の相変化量である。
【０１１３】
また、水蒸気質量分率Ｙｗは、液水ｓと同様に、電界膜内に入り込めないため、電界膜内では水蒸気質量分率Ｙｗ＝０である。
【０１１４】
次に、電界質内の含水量λの輸送方程式について説明する。以下の式（２５）は、電解膜における電界質内の含水量λの輸送方程式である。
・・・（２５）
また、以下の式（２６）は、触媒層内における電解質内の含水量λの輸送方程式である。
・・・（２６）
【０１１５】
また、上記式（２５）および上記式（２６）内の各係数は、以下の式（２７）および式（２８）の関係を有している。
・・・（２７）
・・・（２８）
【０１１６】
ここで、ρｍは電解膜の密度、ＥＷは電解膜の等価質量、Ｖｅｘは電解膜の膨張係数、Ｄｗｌは電解膜の拡散係数、Ｄｉは拡散係数、ｎｄは電気浸透抵抗係数、ｉは電流、Ｆはファラデー定数である。
【０１１７】
なお、電界膜は固体であり液水ｓは入り込めないため、電界膜内では液水ｓ＝０である。
【０１１８】
以上のように、高分子形燃料電池の問題においては、液水ｓは触媒層と多孔質層の二つの領域を移動する物理量であり、各領域にいて異なる輸送方程式により示すことができる。また、水蒸気質量分率Ｙｗは触媒層と多孔質層の二つの領域を移動する物理量であり、各領域において異なる輸送方程式により示すことができる。さらに、電界質内の含水量λは電解質と触媒層の二つの領域を移動する物理量であり、各領域において異なる輸送方程式により示すことができる。
【０１１９】
すなわち、高分子形燃料電池内の液水ｓ、水蒸気質量分率Ｙｗおよび電界質内の含水量λは、それぞれ性質の異なる第一領域および第二領域が第一領域と第二領域との間の境界を介して連結されており、各々の領域を移動する物理量が第一領域および第二領域で異なる輸送方程式によって表されるものであるため、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法によって処理可能な問題である。
【０１２０】
（２）液水ｓの輸送方程式の解析について
まず、液水ｓの輸送方程式の解析について説明する。本実施例１では、第一領域を触媒層、第二領域を多孔質層とした。
【０１２１】
液水ｓの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式（２９）で表される。
・・・（２９）
また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式（３０）で表される。
・・・（３０）
ただし、第一領域側のＰｃ，ＣＬは、以下の式（３１）および式（３２）によって表されるものである。
・・・（３１）
・・・（３２）
【０１２２】
なお、Ｐｃはｓの三次の多項式であり、逆関数Ｐｃ−１を直接求めることが難しい。そこで、以下の方法でＰｃの線形化を行った。
【０１２３】
まず、ｓ＝ａにおいてテイラー展開を行うと、以下の式（３３）で表される一次多項式が得られる。
・・・（３３）
上記式（３３）をＰｃに適用すると、以下の式（３４）が得られる。
・・・（３４）
【０１２４】
以上を踏まえると、上記式（３１）は線形化され、以下の式（３５）として表される。
・・・（３５）
また、第二領域側のＰｃ，ＧＤＬにも適用すると、Ｐｃ，ＧＤＬは以下の式（３６）として表される。
・・・（３６）
ただし、上記式（３５）および上記式（３６）は以下の関係を有する。
【０１２５】
ここで、ａは前時刻のｓ（ａ＝ｓｋ−１）を用いることとすると、上記式（３０）は、以下の式（３７）として表される。
・・・（３７）
【０１２６】
第一領域側の輸送方程式である上記式（１７）を離散化すると、以下の式（３８）として表される。
・・・（３８）
また、第二領域側の輸送方程式である上記式（１８）を離散化すると、以下の式（３９）として表される。
・・・（３９）
なお、上記式（３８）および上記式（３９）は、時間発展計算をオイラー陽解法で解く場合の式である。
【０１２７】
次に、各境界条件についても離散化を行う。以下の式（４０）は、ノイマン条件における上記式（２９）を離散化したものである。
・・・（４０）
また、以下の式（４１）は、ディレクレ条件における上記式（３７）を離散化したものである。
・・・（４１）
【０１２８】
次に、以上の各式を用いて、第一領域および第二領域における液水ｓの輸送方程式の解析について、特に、第一領域内解析手段７および第二領域内解析手段８の解析手順に従って、以下に説明する。
【０１２９】
第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（２９）により境界流束生成量γを算出する。本件の場合は、図１４に示すように、境界における流束（Ｊｆｌｕｘ）は、化学反応等により変化しないと仮定しているため、境界流束生成量γ＝０である。
【０１３０】
次に、第一領域側境界流束算出部７２は、上記式（４０）を用いて第一領域の境界近傍の流束ｆ−（ｓＣＬｋ＋１）を算出する。なお、本件における流束は、各領域における多孔率εに影響されるため、α＝εＣＬであり、β＝εＧＤＬである。
【０１３１】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２より算出された流束ｆ−（ｓＣＬｋ＋１）を用いて上記式（３８）を解き、液水ｓＣＬｋ＋１を求める。
【０１３２】
そして、第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた液水ｓＣＬｋ＋１の値を用いて、第一領域の境界近傍における液水ｓの値ｓ−ｋ＋１を算出する。
【０１３３】
以上により、第一領域内の所定時間Δｔ経過後における液水ｓの値を全て求めることができる。
【０１３４】
次に、第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（３７）により境界における液水ｓの境界変化量ｃを算出する。本件では、境界におけるＰｃの値は化学反応等により変化しないとしているため、境界変化量ｃ＝０である。
【０１３５】
第二領域側境界物質量算出部８２は、上記式（４１）を用いて第二領域の境界近傍の液水ｓＣＬｋ＋１を算出する。なお、本件では、ａ＝１でありｂ＝１である。
【０１３６】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２より算出された液水ｓＣＬｋ＋１を用いて上記式（３９）を解き、液水ＧＤＬｋ＋１を求める。
【０１３７】
そして、第二領域側境界流束算出部は、第二領域内物理量解析部８３により求められた液水ｓＧＤＬｋ＋１の値を用いて、第二領域の境界近傍における液水ｓの流束ｇ＋（ｓＧＤＬｋ＋１）を算出する。
【０１３８】
以上により、第二領域内の所定時間Δｔ経過後における液水ｓの値を全て求めることができる。
【０１３９】
（３）水蒸気質量分率Ｙｗの輸送方程式の解析について
次に、水蒸気質量分率Ｙｗの輸送方程式の解析について説明する。本実施例１では、第一領域を触媒層、第二領域を多孔質層とした。
【０１４０】
水蒸気質量分率Ｙｗの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式（４２）で表される。
・・・（４２）
また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式（４３）で表される。
・・・（４３）
【０１４１】
次に、第一領域側の輸送方程式である上記式（２２）を離散化すると、以下の式（４４）として表される。
・・・（４４）
また、第二領域側の輸送方程式である上記式（２３）を離散化すると、以下の式（４５）として表される。
・・・（４５）
なお、上記式（４４）および上記式（４５）は、時間発展計算をオイラー陽解法で解く場合の式である。
【０１４２】
次に、各境界条件についても離散化を行う。以下の式（４６）は、ノイマン条件における上記式（４２）を離散化したものである。
・・・（４６）
また、以下の式（４７）は、ディレクレ条件における上記式（４３）を離散化したものである。
・・・（４７）
【０１４３】
次に、以上の各式を用いて、第一領域および第二領域における水蒸気質量分率Ｙｗの輸送方程式の解析について、特に、第一領域内解析手段７および第二領域内解析手段８の解析手順に従って、以下に説明する。
【０１４４】
第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（４２）により境界流束生成量γを算出する。本件の場合は、図１４に示すように、境界における流束（Ｊｆｌｕｘ）は、化学反応等により変化しないと仮定しているため、境界流束生成量γ＝０である。
【０１４５】
次に、第一領域側境界流束算出部７２は、上記式（４６）を用いて第一領域の境界近傍の流束ｆ−（ＹＣＬｋ＋１）を算出する。なお、本件における流束は、各領域における多孔率εに影響されるため、α＝εＣＬであり、β＝εＧＤＬである。
【０１４６】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２より算出された流束ｆ−（ｓＣＬｋ＋１）を用いて上記式（４４）を解き、水蒸気質量分率ＹＣＬｋ＋１を求める。
【０１４７】
そして、第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた水蒸気質量分率ＹＣＬｋ＋１の値を用いて、第一領域の境界近傍における水蒸気質量分率Ｙの値Ｙ−ｋ＋１を算出する。
【０１４８】
以上により、第一領域内の所定時間Δｔ経過後における水蒸気質量分率Ｙｗの値を全て求めることができる。
【０１４９】
次に、第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（４３）により境界における水蒸気質量分率Ｙｗの境界変化量ｃを算出する。本件では、境界におけるＹｗの値は化学反応等により変化しないと仮定しているため、境界変化量ｃ＝０である。
【０１５０】
第二領域側境界物質量算出部８２は、上記式（４７）を用いて第二領域の境界近傍の水蒸気質量分率ＹＣＬｋ＋１を算出する。なお、本件では、ａ＝１でありｂ＝１である。
【０１５１】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２より算出された水蒸気質量分率ＹＣＬｋ＋１を用いて上記式（４５）を解き、水蒸気質量分率ＹＧＤＬｋ＋１を求める。
【０１５２】
そして、第二領域側境界流束算出部は、第二領域内物理量解析部８３により求められた水蒸気質量分率ＹＧＤＬｋ＋１の値を用いて、第二領域の境界近傍における水蒸気質量分率Ｙｗの流束ｇ＋（ＹＧＤＬｋ＋１）を算出する。
【０１５３】
以上により、第二領域内の所定時間Δｔ経過後における水蒸気質量分率Ｙｗの値を全て求めることができる。
【０１５４】
（４）電界質内の含水量λの輸送方程式の解析について
次に、電界質内の含水量λの輸送方程式の解析について説明する。本実施例１では、第一領域を電解膜、第二領域を触媒層とした。
【０１５５】
電界質内の含水量λの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式（４８）で表される。
・・・（４８）
また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式（４９）で表される。
・・・（４９）
【０１５６】
次に、第一領域側の輸送方程式である上記式（２５）を離散化すると、以下の式（５０）として表される。
・・・（５０）
また、第二領域側の輸送方程式である上記式（２６）を離散化すると、以下の式（５１）として表される。
・・・（５１）
なお、上記式（５０）および上記式（５１）は、時間発展計算をオイラー陽解法で解く場合の式である。
【０１５７】
次に、各境界条件についても離散化を行う。以下の式（５２）は、ノイマン条件における上記式（４８）を離散化したものである。
・・・（５２）
また、以下の式（５３）は、ディレクレ条件における上記式（４９）を離散化したものである。
・・・（５３）
【０１５８】
次に、以上の各式を用いて、第一領域および第二領域における電界質内の含水量λの輸送方程式の解析について、特に、第一領域内解析手段７および第二領域内解析手段８の解析手順に従って、以下に説明する。
【０１５９】
第一領域内解析手段７では、まず、境界流束生成量算出部７１が、上記式（４８）により境界流束生成量γを算出する。本件の場合は、図１４に示すように、境界における流束（Ｊｆｌｕｘ）は、化学反応等により変化しないと仮定しているため、境界流束生成量γ＝０である。
【０１６０】
次に、第一領域側境界流束算出部７２は、上記式（４６）を用いて第一領域の境界近傍の流束ｆ−（λＣＬｋ＋１）を算出する。なお、本件では、α＝１であり、β＝１である。
【０１６１】
第一領域内物理量解析部７３は、第一領域側境界流束算出部７２より算出された流束ｆ−（λＣＬｋ＋１）を用いて上記式（５２）を解き、電界質内の含水量λＣＬｋ＋１を求める。
【０１６２】
そして、第一領域側境界物理量算出部７４は、第一領域内物理量解析部７３により求められた電界質内の含水量λＣＬｋ＋１の値を用いて、第一領域の境界近傍における電界質内の含水量λの値λ−ｋ＋１を算出する。
【０１６３】
以上により、第一領域内の所定時間Δｔ経過後における電界質内の含水量λの値を全て求めることができる。
【０１６４】
次に、第二領域内解析手段８では、まず、境界変化量算出部８１が、上記式（４９）により境界における電界質内の含水量λの境界変化量ｃを算出する。本件では、境界におけるλの値は化学反応等により変化しないと仮定しているため、境界変化量ｃ＝０である。
【０１６５】
第二領域側境界物質量算出部８２は、上記式（５３）を用いて第二領域の境界近傍の電界質内の含水量λＣＬｋ＋１を算出する。なお、本件では、ａ＝１でありｂ＝１である。
【０１６６】
第二領域内物理量解析部８３は、第二領域側境界物理量算出部８２より算出された電界質内の含水量λＣＬｋ＋１を用いて上記式（５１）を解き、電界質内の含水量λＧＤＬｋ＋１を求める。
【０１６７】
そして、第二領域側境界流束算出部８４は、第二領域内物理量解析部８３により求められた電界質内の含水量λＧＤＬｋ＋１の値を用いて、第二領域の境界近傍における電界質内の含水量λの流束ｇ＋（λＧＤＬｋ＋１）を算出する。
【０１６８】
以上により、第二領域内の所定時間Δｔ経過後における電界質内の含水量λｗの値を全て求めることができる。
【実施例２】
【０１６９】
『本発明における薄膜近似を用いた高分子形燃料電池の解析』
本実施例２では、上述した第二実施形態に基づく、本発明における薄膜近似を用いた高分子形燃料電池の解析手法について説明する。
【０１７０】
本実施例２における解析場は、図１５に示すように、電解膜、触媒層および多孔質層の３つの領域からなり、中間に位置する触媒層における物理量を薄膜近似して処理するものである。
【０１７１】
（１）各層における輸送方程式について
各層における輸送方程式について、図１５を参照しながら説明する。図１５に示すように、電解膜および多孔質層における液水ｓ、水蒸気質量分率Ｙｗおよび電解質内の含水量λの輸送方程式は、上述した実施例１におけるものと同一である。
【０１７２】
具体的には、電解膜内では、液水ｓ＝０、水蒸気質量分率Ｙｗ＝０であるとともに、電解質内の含水量λの輸送方程式は上記式（２５）により表される。また、多孔質層内では、液水ｓの輸送方程式は上記式（１８）により表され、水蒸気質量分率Ｙｗの輸送方程式は上記式（２３）により表されるとともに、電解質内の含水量λ＝０である。
【０１７３】
なお、以下の説明においては、上記式（２５）は、基礎方程式形で表された、以下の式（５４）として表した。
・・・（５４）
同様に、上記式（１８）は、以下の式（５５）とした。
・・・（５５）
さらに、上記式（２３）は、以下の式（５６）とした。
・・・（５６）
【０１７４】
また、本実施例２では、触媒層内の物理量の輸送は一定であると仮定して、その平均値ｓ、Ｙｗ、λおよび平均流束ｈ（ｓ）、ｈ（Ｙｗ）、ｈ（λ）として定義した。
【０１７５】
（２）各輸送方程式の解析について
本実施例２では、液水ｓ、水蒸気質量分率Ｙｗおよび電解質内の含水量λを同時に処理する。また、多孔質層を第一領域、触媒層を第二領域および電解膜を第三領域とした。
【０１７６】
電解膜から多孔質層への物理量の移動に対しては、ノイマン条件を与えた。ただし、電解膜における電解質内の含水量λは、実際には多孔質層内の水蒸気質量分率Ｙｗおよび液水ｓとして移動するが、すべてが多孔質層内の水蒸気質量分率Ｙｗとして移動するものとして仮定した。なお、この仮定下においては、多孔質層内の水蒸気質量分率Ｙｗが飽和量を超えた際には直ちに液水ｓに変換される関係にあり、電解質内の含水量λは間接的に液水ｓになるものであるため、物理的な関連性は保たれているものと考えられる。
【０１７７】
以上の仮定を踏まえると、第三領域側境界におけるノイマン条件は、以下の式（５７）
で表される。・・・（５７）
また、第一領域側境界におけるノイマン条件は、以下の式（５８）で表される。
・・・（５８）
【０１７８】
また、多孔質層から電解膜への物理量の移動に対しては、ディレクレ条件を与えた。ここで、触媒層内における電解質内の含水量λは、平衡状態に達していると考え、以下の式（５９）に示すように、液水ｓとａｗ（water activity）の関数で表されるものとした。
・・・（５９）
【０１７９】
第三領域側境界におけるディレクレ条件は、以下の式（６０）で表される。
・・・（６０）
また、第一領域境界においては、液水ｓおよび水蒸気質量分率Ｙｗのディレクレ条件が与えられ、液水ｓに関しては、以下の式（６１）で表される。
・・・（６１）
また、水蒸気質量分率Ｙｗに関しては、以下の式（６２）で表される。
・・・（６２）
【０１８０】
次に、上述した各式の離散化を行う。第一領域における電解質内の含水量λの輸送方程式である上記式（５４）の離散化をすると、以下の式（６３）として表される。
・・・（６３）
また、第三領域における液水ｓの輸送方程式である上記式（５５）の離散化をすると、以下の式（６４）として表される。
・・・（６４）
さらに、第三領域における水蒸気質量分率Ｙｗの輸送方程式である上記式（５６）を離散化すると、以下の式（６５）として表される。
・・・（６５）
【０１８１】
次に、各境界条件についても離散化を行う。第三領域側境界におけるノイマン条件の上記式（５７）を離散化すると、以下の式（６６）として表される。
・・・（６６）
また、本実施例２では、電解質内の含水量λは、すべてが多孔質層内の水蒸気質量分率Ｙｗとして移動するものと仮定したため、電解質内の含水量λと水蒸気質量分率Ｙｗとの間には、以下の式（６７）の関係が成り立つ。
・・・（６７）
また、第一領域側境界におけるノイマン条件の上記式（５８）の離散化を行うと、以下の式（６８）として表される。
・・・（６８）
【０１８２】
次に、第三領域側境界のディレクレ条件の上記式（６０）を離散化すると、以下の式（６９）として表される。
・・・（６９）
また、第三領域側境界における液水ｓのディレクレ条件の上記式（６１）の離散化をすると、以下の式（７０）として表される。
・・・（７０）
さらに、第一領域側境界における水蒸気質量分率Ｙｗのディレクレ条件の上記式（６２）を離散化すると、以下の式（７１）として表される。
・・・（７１）
【０１８３】
また、電解質内の含水量λを表した上記式（５９）を離散化すると、以下の式（７２）として表される。
・・・（７２）
【０１８４】
次に、以上の各式を用いて、第一領域および第三領域における各物理量の輸送方程式の解析について、特に、第二実施形態における第一領域内解析手段９および第三領域内解析手段１０の解析手順に従って、以下に説明する。
【０１８５】
第一領域内解析手段９では、まず、第一領域側境界流束生成量算出部９１が、境界流束生成量γ＋を算出する。本件では、境界流束生成量γ＋は、以下の式（７３）として算出される。
・・・（７３）
なお、陰極側の計算には、以下の式（７４）で表されるバトラー・ボルマー方程式を用いた。
・・・（７４）
ここで、電子２ｍｏｌ流れたときにＨ２Ｏは１ｍｏｌ発生するため、以下の式（７５）で表される水が発生する。
・・・（７５）
この式（７５）に、触媒層の厚みＬを掛けると、上記式（７３）の陰極側の値が導き出される。
【０１８６】
次に、第三領域側境界境界流束生成量算出部９２が、境界流束生成量γ−を算出する。本件では、境界流束生成量γ−＝０である。
【０１８７】
次に、第一領域側境界流束算出部９３は、上記式（６６）、上記式（６７）および上記式（６８）を用いて第一領域側境界の流束ｇ（Ｙｗ、ＧＤＬｋ＋１）を算出する。なお、本件におけるαおよびβを以下に示す。
【０１８８】
第一領域内物理量解析部９４は、第一領域側境界流束算出部７３より算出された流束ｇ（Ｙｗ、ＧＤＬｋ＋１）を用いて上記式（６５）を解き、多孔質層内の水蒸気質量分率Ｙｊｋ＋１を求める。
【０１８９】
そして、第一領域側境界物理量算出部９５は、第一領域内物理量解析部９４により求められた多孔質層内の水蒸気質量分率Ｙｊｋ＋１の値を用いて、第一領域側境界における水蒸気質量分率Ｙｗの値Ｙ＋ｋ＋１を算出する。
【０１９０】
以上により、第一領域内の所定時間Δｔ経過後における水蒸気質量分率Ｙｗの値を全て求めることができる。
【０１９１】
次に、第三領域内解析手段１０では、第三領域側境界変化量算出部１０１が、上記式（６０）により境界変化量ｃ−を算出する。本件では、境界変化量ｃ−＝０である。
【０１９２】
また、第一領域側境界変化量算出部１０２が、上記式（６１）および上記式（６２）により、境界変化量ｃｓ，＋および境界変化量ｃＹ，＋を算出する。本件では、境界変化量ｃｓ，＋＝０であり、境界変化量ｃＹ，＋＝０である。
【０１９３】
次に、第三領域側境界物質量算出部１０３は、上記式（６９）、上記式（７１）および上記式（７２）を用いて第三領域側境界の電界質内の含水量λ−ｋ＋１を算出する。なお、本件では、ａおよびｂに関する値はすべて１である。
【０１９４】
第三領域内物理量解析部１０４は、第三領域側境界物質量算出部１０３より算出された電界質内の含水量λ−ｋ＋１を用いて上記式（６３）を解き、電界質内の含水量λｉｋ＋１を求める。
【０１９５】
そして、第三領域側境界流束算出部１０５は、第三領域内物理量解析部１０４により求められた電界質内の含水量λｉｋ＋１の値を用いて、第三領域側境界における電界質内の含水量λの流束ｆ＋（λＧＤＬｋ＋１）を算出する。
【０１９６】
以上により、第三領域内の所定時間Δｔ経過後における電界質内の含水量λｗの値を全て求めることができる。
【０１９７】
なお、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法による解析対象は、前述した高分子燃料電池の問題に限定されるものではなく、他にも適用することができる。例えば、二次電池におけるイオン伝導問題や多孔質内における燃焼問題等に使用することができる。
【符号の説明】
【０１９８】
１ コンピュータ
１ａ，１ｂ 輸送方程式解析プログラム
２ 入力手段
３ 表示手段
４ 記憶手段
５ 演算手段
６ 物理量取得手段
７ 第一領域内解析手段
８ 第二領域内解析手段
９ 第一領域内解析手段
１０ 第三領域内解析手段
４１ プログラム記憶部
４２ 物理量記憶部
７１ 境界流束生成量算出部
７２ 第一領域側境界流束算出部
７３ 第一領域内物理量解析部
７４ 第一領域側境界物理量算出部
８１ 境界変化量算出部
８２ 第二領域側境界物理量算出部
８３ 第二領域内物理量解析部
８４ 第二領域側境界流束算出部
９１ 第一領域側境界流束生成量算出部
９２ 第三領域側境界流束生成量算出部
９３ 第一領域側境界流束算出部
９４ 第一領域内物理量解析部
９５ 第一領域側境界物理量算出部
１０１ 第三領域側境界変化量算出部
１０２ 第一領域側境界変化量算出部
１０３ 第三領域側境界物理量算出部
１０４ 第三領域内物理量解析部
１０５ 第三領域側境界流束算出部
【特許請求の範囲】
【請求項１】
それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析プログラムであって、
前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、
前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、
前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析手段と
してコンピュータを機能させるとともに、
前記第一領域内解析手段は、
前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出部と、
前記境界流束生成量算出部により算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出部と、
前記第一領域側境界流束算出部により算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、
前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出部として機能させ、
前記第二領域内解析手段は、
前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出部と、
前記境界変化量算出部により算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出部と、
前記第二領域側境界物理量算出部により算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析部と、
前記第二領域内物理量解析部により求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出部
として機能させる輸送方程式解析プログラム。
【請求項２】
それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析プログラムであって、
前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、
前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、
前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析手段と
してコンピュータを機能させるとともに、
前記第一領域内解析手段は、
前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出部と、
前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出部と、
前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出部と、
前記第一領域側境界流束算出部により算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、
前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出部として機能させ、
前記第三領域内解析手段は、
前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出部と、
前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出部と、
前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出部と、
前記第三領域側境界物理量算出部により算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析部と、
前記第三領域内物理量解析部により求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出部
として機能させる輸送方程式解析プログラム。
【請求項３】
それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析方法であって、
前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得ステップと、
前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析ステップと、
前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析ステップと
を有しており、
前記第一領域内解析ステップは、
前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出ステップと、
前記境界流束生成量算出ステップにより算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出ステップと、
前記第一領域側境界流束算出ステップにより算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析ステップと、
前記第一領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出ステップと
を有するとともに、
前記第二領域内解析ステップは、
前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出ステップと、
前記境界変化量算出ステップにより算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出ステップと、
前記第二領域側境界物理量算出ステップにより算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析ステップと、
前記第二領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出ステップ
とを有する輸送方程式解析方法。
【請求項４】
それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析方法であって、
前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得ステップと、
前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析ステップと、
前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析ステップと
を有しており、
前記第一領域内解析ステップは、
前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出ステップと、
前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出ステップと、
前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出ステップと、
前記第一領域側境界流束算出ステップにより算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析ステップと、
前記第一領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出ステップと
を有するとともに、
前記第三領域内解析ステップは、
前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出ステップと、
前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出ステップと、
前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出ステップと、
前記第三領域側境界物理量算出ステップにより算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析ステップと、
前記第三領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出ステップと
を有する輸送方程式解析方法。
【請求項１】
それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析プログラムであって、
前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、
前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、
前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析手段と
してコンピュータを機能させるとともに、
前記第一領域内解析手段は、
前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出部と、
前記境界流束生成量算出部により算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出部と、
前記第一領域側境界流束算出部により算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、
前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出部として機能させ、
前記第二領域内解析手段は、
前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出部と、
前記境界変化量算出部により算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出部と、
前記第二領域側境界物理量算出部により算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析部と、
前記第二領域内物理量解析部により求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出部
として機能させる輸送方程式解析プログラム。
【請求項２】
それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析プログラムであって、
前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、
前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、
前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析手段と
してコンピュータを機能させるとともに、
前記第一領域内解析手段は、
前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出部と、
前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出部と、
前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出部と、
前記第一領域側境界流束算出部により算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、
前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出部として機能させ、
前記第三領域内解析手段は、
前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出部と、
前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出部と、
前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出部と、
前記第三領域側境界物理量算出部により算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析部と、
前記第三領域内物理量解析部により求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出部
として機能させる輸送方程式解析プログラム。
【請求項３】
それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析方法であって、
前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得ステップと、
前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析ステップと、
前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析ステップと
を有しており、
前記第一領域内解析ステップは、
前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出ステップと、
前記境界流束生成量算出ステップにより算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出ステップと、
前記第一領域側境界流束算出ステップにより算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析ステップと、
前記第一領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出ステップと
を有するとともに、
前記第二領域内解析ステップは、
前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出ステップと、
前記境界変化量算出ステップにより算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出ステップと、
前記第二領域側境界物理量算出ステップにより算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析ステップと、
前記第二領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出ステップ
とを有する輸送方程式解析方法。
【請求項４】
それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析方法であって、
前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得ステップと、
前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析ステップと、
前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析ステップと
を有しており、
前記第一領域内解析ステップは、
前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出ステップと、
前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出ステップと、
前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出ステップと、
前記第一領域側境界流束算出ステップにより算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析ステップと、
前記第一領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出ステップと
を有するとともに、
前記第三領域内解析ステップは、
前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出ステップと、
前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出ステップと、
前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出ステップと、
前記第三領域側境界物理量算出ステップにより算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析ステップと、
前記第三領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出ステップと
を有する輸送方程式解析方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【公開番号】特開２０１３−６１６９１（Ｐ２０１３−６１６９１Ａ）
【公開日】平成２５年４月４日（２０１３．４．４）
【国際特許分類】
【出願番号】特願２０１１−１９７７５５（Ｐ２０１１−１９７７５５）
【出願日】平成２３年９月１０日（２０１１．９．１０）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成１９年度〜平成２０年度 独立行政法人新エネルギー・産業技術総合開発機構 固体高分子形燃料電池実用化戦略的技術開発／次世代燃料電池技術開発委託研究、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（５０４１７３４７１）国立大学法人北海道大学 (971)
【公開日】平成２５年４月４日（２０１３．４．４）
【国際特許分類】
【出願日】平成２３年９月１０日（２０１１．９．１０）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成１９年度〜平成２０年度 独立行政法人新エネルギー・産業技術総合開発機構 固体高分子形燃料電池実用化戦略的技術開発／次世代燃料電池技術開発委託研究、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（５０４１７３４７１）国立大学法人北海道大学 (971)
[ Back to top ]