"それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析プログラムであって、前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出部と、前記境界流束生成量算出部により算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第二領域内解析手段は、前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出部と、前記境界変化量算出部により算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出部と、前記第二領域側境界物理量算出部により算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析部と、前記第二領域内物理量解析部により求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出部として機能させる輸送方程式解析プログラム。それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析プログラムであって、前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出部と、前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出部と、前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第三領域内解析手段は、前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出部と、前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出部と、前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出部と、前記第三領域側境界物理量算出部により算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析部と、前記第三領域内物理量解析部により求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出部として機能させる輸送方程式解析プログラム。それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析方法であって、前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得ステップと、前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析ステップと、前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析ステップとを有しており、前記第一領域内解析ステップは、前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出ステップと、前記境界流束生成量算出ステップにより算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出ステップと、前記第一領域側境界流束算出ステップにより算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析ステップと、前記第一領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出ステップとを有するとともに、前記第二領域内解析ステップは、前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出ステップと、前記境界変化量算出ステップにより算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出ステップと、前記第二領域側境界物理量算出ステップにより算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析ステップと、前記第二領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出ステップとを有する輸送方程式解析方法。それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析方法であって、前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得ステップと、前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析ステップと、前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析ステップとを有しており、前記第一領域内解析ステップは、前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出ステップと、前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出ステップと、前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出ステップと、前記第一領域側境界流束算出ステップにより算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析ステップと、前記第一領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出ステップとを有するとともに、前記第三領域内解析ステップは、前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出ステップと、前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出ステップと、前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出ステップと、前記第三領域側境界物理量算出ステップにより算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析ステップと、前記第三領域内物理量解析ステップにより求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出ステップとを有する輸送方程式解析方法。しかしながら、特許文献1に記載された発明を含む、これまでの燃料電池の数値解析プログラムは以下のような問題点を有している。すなわち、1.従来の数値解析プログラムは、発電セルの平面内における所定の物理量の分布を解析するものである。よって、それぞれ性質の異なる多層領域間の物理量の移動を正確に求めているわけではなく、多層領域間を移動する物理量の計算や当該物理量の計算に基づいた発電量の計算等の解析結果がその精度において不十分であった。2.燃料電池などの多層領域間における物理量は、各層の境界における化学反応等により各種の物理量の変化を伴い、解析結果は前記境界における境界条件に強い影響を受ける。そのため、与えた境界条件や初期値によっては収束性が悪くなったり、解が得られない場合があった。3.燃料電池の解析において、多孔質材と電解質膜との間にある触媒層は、他の層に比べて非常に薄く、解析精度は触媒層近傍における格子解像度に依存していた。そのため、格子を十分に小さくする必要があったが、その分、計算負荷が大きくなっていた。本発明は、このような問題点を解決するためになされたものであって、解析精度の向上、解析速度の向上、解析負荷の軽減ならびに収束性の安定を図ることのできる輸送方程式解析プログラムおよび輸送方程式解析方法を提供することを目的としている。本発明に係る輸送方程式解析方法を用いた輸送方程式解析プログラムは、それぞれ性質の異なる第一領域および第二領域が前記第一領域と前記第二領域との間の境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域および前記第二領域で異なる輸送方程式によって表される場合の前記各輸送方程式の数値解を求めるための輸送方程式解析プログラムであって、前記第一領域および前記第二領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記境界へ流入出する前記物理量の流束と前記第二領域側から前記境界へ流入出する前記物理量の流束との連続性を示す流束境界条件式、および前記第一領域の輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す物理量境界条件式、および前記第二領域の輸送方程式に基づいて前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界流束生成量を算出する境界流束生成量算出部と、前記境界流束生成量算出部により算出された境界流束生成量および前記第二領域における前記境界近傍の前記流束を前記流束境界条件式に代入して前記所定時間後の前記第一領域における前記境界近傍の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記境界近傍の流束を前記第一領域の前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域における前記境界近傍の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第二領域内解析手段は、前記物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記境界における前記物理量の境界変化量を算出する境界変化量算出部と、前記境界変化量算出部により算出された境界変化量および前記第一領域における前記境界近傍の前記物理量を前記物理量境界条件式に代入して前記所定時間後の前記第二領域における前記境界近傍の物理量を算出する第二領域側境界物理量算出部と、前記第二領域側境界物理量算出部により算出された前記境界近傍の物理量を前記第二領域の前記輸送方程式に代入して前記所定時間後の前記第二領域内における前記物理量の数値解を求める第二領域内物理量解析部と、前記第二領域内物理量解析部により求められた数値解から前記所定時間後の前記第二領域における前記境界近傍の流束を算出する第二領域側境界流束算出部として機能させる。すなわち、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法は、第一領域と第二領域との境界において、ディレクレ条件に基づく境界条件とノイマン条件に基づく境界条件の2つの境界条件式を用いるとともに、その2つの境界条件式を所定の順序で解くことによって、汎用的な輸送方程式解析プログラムおよび輸送方程式解析方法を構築したものである。また、本発明に係る輸送方程式解析方法を用いた輸送方程式解析プログラムは、それぞれ性質の異なる第一領域、第二領域および第三領域において、前記第一領域と前記第二領域とがこれらの間の第一領域側境界を介して連結されているとともに、前記第二領域と前記第三領域とがこれらの間の第三領域側境界を介して連結されており、前記各々の領域を移動する所定の物理量が前記第一領域と前記第三領域とで異なる輸送方程式によって表されるとともに、前記第二領域が薄膜近似とみなせる場合における前記第一領域および前記第三領域の前記各輸送方程式の数値解を求める輸送方程式解析プログラムであって、前記第一領域および第三領域を空間的に離散化した各格子点における前記物理量を取得する物理量取得手段と、前記第一領域側から前記第一領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第一領域側境界へ流入出する前記物理量の流束との連続性を示す第一領域側流束境界条件式、前記第三領域側から前記第三領域側境界へ流入出する前記物理量の流束と前記第二領域側から前記第三領域側境界へ流入出する前記物理量の流束との連続性を示す第三領域側流束境界条件式、および前記第一領域における前記輸送方程式に基づいて所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内解析手段と、前記第一領域側境界における前記第一領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第一領域側物理量境界条件式、前記第三領域側境界における前記第三領域側の前記物理量と前記第二領域側の前記物理量との連続性を示す第三領域側物理量境界条件式、および前記第三領域における前記輸送方程式に基づいて前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内解析手段としてコンピュータを機能させるとともに、前記第一領域内解析手段は、前記第一領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界流束生成量を算出する第一領域側境界流束生成量算出部と、前記第三領域側流束境界条件式において定義された境界流束生成量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界流束生成量を算出する第三領域側境界流束生成量算出部と、前記第三領域側境界における前記境界流束生成量および前記第三領域側境界の前記流束を前記第三領域側流束境界条件式に代入することにより所定時間後の前記第二領域内の平均流束を算出するとともに、この第二領域内の前記平均流束と前記第一領域側境界における前記境界流束生成量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第一領域側境界の流束を算出する第一領域側境界流束算出部と、前記第一領域側境界流束算出部により算出された前記第一領域側境界の流束を前記第一領域における前記輸送方程式に代入して前記所定時間後の前記第一領域内における前記物理量の数値解を求める第一領域内物理量解析部と、前記第一領域内物理量解析部により求められた数値解から前記所定時間後の前記第一領域側境界の物理量を算出する第一領域側境界物理量算出部として機能させ、前記第三領域内解析手段は、前記第三領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第三領域側境界における前記物理量の境界変化量を算出する第三領域側境界変化量算出部と、前記第一領域側物理量境界条件式において定義された境界変化量に係る関数に取得した前記物理量を代入して前記第一領域側境界における前記物理量の境界変化量を算出する第一領域側境界変化量算出部と、前記第一領域側境界における前記境界変化量および前記第一領域側境界の前記物理量を前記第一領域側流束境界条件式に代入することにより前記所定時間後の前記第二領域内における前記物理量の平均値を算出するとともに、この第二領域内の前記平均値と前記第三領域側境界における前記物理量の境界変化量を前記第三領域側流束境界条件式に代入することにより前記所定時間後の前記第三領域側境界の前記物理量を算出する第三領域側境界物理量算出部と、前記第三領域側境界物理量算出部により算出された前記第三領域側境界の物理量を前記第三領域側の前記輸送方程式に代入して前記所定時間後の前記第三領域内における前記物理量の数値解を求める第三領域内物理量解析部と、前記第三領域内物理量解析部により求められた数値解から前記所定時間後の前記第三領域側境界の前記流束を算出する第三領域側境界流束算出部として機能させる。すなわち、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法は、第二領域が薄膜近似と捉えられる場合であって、前記第一領域側境界および前記第三領域側境界に、ディレクレ条件に基づく境界条件とノイマン条件に基づく境界条件の2つの境界条件式を用いるとともに、それらの合計4つの境界条件式を所定の順序で解くことによって、薄膜近似を用いた汎用的な輸送方程式解析プログラムおよび輸送方程式解析方法を構築したものである。本発明によれば、解析精度の向上、解析速度の向上、解析負荷の軽減ならびに収束性の安定を図ることができる。以下、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法の第一実施形態について図面を用いて説明する。図1は、本第一実施形態における輸送方程式解析プログラム1aを備えたコンピュータの構成を示すブロック図である。また、図2は、本第一実施形態で解析を行う第一領域および第二領域を示す概略図である。前記第一領域および前記第二領域は、境界を介して連結されており、それぞれ性質の異なる構造および/または物質からなる。また、前記各領域の間を移動する所定の物理量は、異なる輸送方程式により表される。例えば、第一領域における輸送方程式は、以下の式1により表される。式1また、第二領域における輸送方程式は、以下の式2により表される。式2ここで、上記式1におけるφは第一領域における物理量であり、上記式2におけるψは第二領域における物理量である。なお、本発明において、上記φと上記ψは同じ物理量を表すものであり、上記式1および上記式2における物理量を同じ符号を用いて記載することも可能であるが、以下に記載される各式における物理的意味を明確にするため、第一領域における物理量と第二領域における物理量を別表記とした。また、本発明における物理量は、例えば、水、水蒸気および熱等を示すものであり、各領域において輸送方程式により表されるもの全ての物理量である。本第一実施形態におけるコンピュータ1Aは、図1に示すように、主として、入力手段2、表示手段3、記憶手段4および演算処理手段5から構成されている。以下、各構成について詳細に説明する。入力手段2は、テキストや数値、記号等を入力する操作キー、操作マウス等からなる。本第一実施形態においては、輸送方程式解析プログラム1aを実行する際に用いられる初期条件や設定条件等の数値の入力操作等に用いることができるようになっている。表示手段3は、画像やテキストデータを表示する液晶ディスプレーやCRTディスプレー、タッチパネル等からなり、輸送方程式解析プログラム1aにおけるユーザーインターフェースとして、入力手段2により入力された内容や解析結果等を表示できるようになっている。記憶手段4は、ROMReadOnlyMemory、RAMRandomAccessMemory、ハードディスク、フラッシュメモリ等によって構成されており、各種のデータを記憶するとともに、演算処理手段5が演算を行う際のワーキングエリアとして機能するものである。本第一実施形態において、記憶手段4は、図1に示すように、主として、輸送方程式解析プログラム1aを記憶するプログラム記憶部41と、空間的に離散化された各格子点における物理量を記憶させておく物理量記憶部42とを有する。なお、本第一実施形態における物理量記憶部42には、各格子点における物理量とともに、各格子点における流束が記憶されている。演算処理手段5は、CPUCentralProcessingUnit等から構成されており、記憶手段4にインストールされた第一実施形態の輸送方程式解析プログラム1aを実行させることにより、図1に示すように、物理量取得手段6、第一領域内解析手段7および第二領域内解析手段8としてコンピュータ1Aを機能させるようになっている。以下、各構成手段についてより詳細に説明する。物理量取得手段6は、記憶手段4の物理量記憶部42に記憶された所定の時間tの各格子点における物理量を取得するものである。また、本第一実施形態における物理量取得手段6は、各格子点における物理量とともに、各格子点の流束も取得している。なお、流束は取得した物理量から計算して求めることも可能である。第一領域内解析手段7は、物理量取得手段6により取得された物理量の時間tにより、所定時間Δt後における第一領域内の物理量の数値解を求めるものである。例えば、物理量取得手段6により取得された物理量が初期値である場合、物理量取得手段6により取得された物理量の時間tは、t=0であり、第一領域内解析手段7は、その時間t=0から所定時間Δt後の物理量を求めるものである。本第一実施形態における第一領域内解析手段7は、図1に示すように、境界流束生成量算出部71と、第一領域側境界流束算出部72と、第一領域内物理量解析部73と、第一領域側境界物理量算出部74とを有している。境界流束生成量算出部71は、以下の式3の流束境界条件式により定義された境界流束生成量γに係る関数γφ,ψを用いて前記境界流束生成量γの値を算出するものである。3よって、境界流束生成量算出部71では、前記境界流束生成量γの関数γφ,ψに、前記物理量取得手段6により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界流束生成量γを算出している。なお、この境界流束生成量γは、境界における化学反応等により物理量が変化したことにより生じる流束の差を示す値であり、算出された境界流束生成量γの符号が正の場合には単位時間当たりの生成量を、符号が負の場合には単位時間当たりの消滅量を表す。また、上記式3で示される流束境界条件式は、第一領域と第二領域との間のノイマン条件による境界条件を表す式であり、前記第一領域側から前記境界へ流入出する前記物理量の流束fφと前記第二領域側から前記境界へ流入出する前記物理量の流束gψとの連続性を表したものである。ここで、αは第一領域における係数であり、βは第二領域における係数である。例えば、物理量が液水であって第一領域が多孔質層である場合、液水は多孔質層内の孔の中しか通過することができない。よって、第一領域における流束は、流束nfφと第一領域の多孔率孔の体積/全体の体積εとの積で表される。よって、この場合はα=εとなる。第一領域側境界流束算出部72は、前記境界流束生成量算出部71により算出された境界流束生成量γ、および前記物理量取得手段により取得された第二領域における境界近傍の流束gψkを、流束境界条件式に代入して所定時間Δt経過後の第一領域における境界近傍の流束f-φk1を算出するものである。よって、第一領域側境界流束算出部72における流束境界条件式は、所定時間Δt経過後の流束f-φk1を算出するため、以下の式4に示すように時間方向に離散化されている。式4ここで、kは時間を表すものであって物理量取得手段6により取得された物理量の有する時間tに対応するものである。また、k1は、時間tから所定時間Δt経過後の時間を表すものである。また、は第二領域における境界近傍の値であることを示しており、具体的には、境界から第二領域側にある一つ目の格子点における値を示している。また、-は第一領域における境界近傍の値であることを示しており、具体的には、境界から第一領域側にある一つ目の格子点における値を示している。よって、式4は、第二領域から境界へと流入する時間kにおける流束gψkが、境界において境界流束生成量γ分だけ増減され、所定時間Δt経過後の時間k1には、第一領域における境界近傍へとf-φk1として流出していることを示している。第一領域内物理量解析部73は、第一領域側境界流束算出部72により算出された境界近傍の流束αnf-φk1を第一領域の輸送方程式に代入して所定時間Δt経過後の第一領域内における物理量φik1の数値解を求めるものである。よって、第一領域内物理量解析部73における第一領域の輸送方程式は、所定時間Δt経過後の第一領域内における物理量の数値解を求めるために空間的に離散化すると、上記式1を以下の式5として表される。5ここで、iは、第一領域内の格子点の位置を表すものであり、本第一実施形態においては、境界から2番目の格子点から、i=1,2,3,と番号が割り当てられる。本第一実施形態における第一領域内物理量解析部73は、第一領域側境界流束算出部72により算出した第一領域における境界近傍の流束f-φk1を上記式5に代入して導かれる連立方程式を解くことで、第一領域内の各格子点における所定時間Δt経過後の物理量φik1を算出する。なお、上記式5による連立方程式を解く方法は特に限定されるものではなく、陰解法や陽解法等から適宜選択されるものである。第一領域側境界物理量算出部74は、第一領域内物理量解析部73により求められた数値解から所定時間Δt経過後の第一領域における境界近傍の物理量aφ-k1を算出するものである。次に、第二領域内解析手段8について説明する。第二領域内解析手段8は、所定時間Δt経過後における第二領域内の物理量の数値解を求めるものである。本第一実施形態における第二領域内解析手段8は、図1に示すように、境界変化量算出部81と、第二領域側境界物理量算出部82と、第二領域内物理量解析部83と、第二領域側境界流束算出部84とを有している。境界変化量算出部81は、以下の式6の物理量境界条件式により定義された境界変化量cに係る関数cφ,ψを用いて前記境界変化量cの値を算出するものである。6よって、境界変化量算出部81では、前記境界変化量cの関数cφ,ψに、前記物理量取得手段6により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界変化量cを算出している。なお、この境界変化量cは、境界における流束の伝達率により生じる第一領域側の物理量と第二領域側の物理量との差を示す値である。また、上記式6で示される物理量境界条件式は、第一領域と第二領域との間のディレクレ条件による境界条件を表す式であり、境界における第一領域側の物理量φと第二領域側の物理量ψとの連続性を示している。第二領域側境界物理量算出部82と境界変化量算出部81とにより算出された境界変化量cおよび第一領域における境界近傍の物理量φ-kを物理量境界条件式に代入して、所定時間Δt経過後の第二領域における境界近傍の物理量ψk1を算出するものである。よって、第二領域側境界物理量算出部82における物理量境界条件式は、所定時間Δt経過後の物理量bψkψk1を算出するため、以下の式7に示すように時間方向に離散化されている。7ここで、aは第一領域における係数であり第一領域の物理量φの関数として与えられている。また、bは第二領域における係数であり第二領域の物理量ψの関数として与えられる。第二領域内物理量解析部83は、第二領域側境界物理量算出部82により算出された境界近傍の物理量bψkψk1を第二領域の輸送方程式に代入して所定時間Δt経過後の第二領域内における物理量ψjk1の数値解を求めるものである。よって、第二領域内物理量解析部83における第二領域の輸送方程式は、所定時間Δt経過後の第二領域内における物理量ψjk1の数値解を求めるために空間的に離散化すると、上記式2を以下の式8として表される。8ここで、jは、第二領域内の格子点の位置を表すものであり、本第一実施形態においては、境界から2番目の格子点から、j=1,2,3,と番号が割り当てられる。本第一実施形態における第二領域内物理量解析部83は、第二領域側境界物理量算出部82により算出した第二領域における境界近傍の物理量ψk1を上記式8に代入して導かれる連立方程式を解くことで、第二領域内の各格子点における所定時間Δt経過後の物理量ψjk1を算出している。なお、上記式8による連立方程式を解く方法は特に限定されるものではなく、陰解法や陽解法等から適宜選択されるものである。第二領域側境界流束算出部84は、第二領域内物理量解析部83により求められた数値解から所定時間Δt経過後の第一領域における境界近傍の流束gψk1を算出するものである。次に、本第一実施形態の輸送方程式解析プログラム1aにおける各構成の作用および輸送方程式解析方法について、図3ないし図5に示すフローチャート図および図6示す各構成の機能を示す機能ブロック図を用いて説明する。図3および図6に示すように、本第一実施形態の輸送方程式解析プログラム1aおよび輸送方程式解析方法では、まず、物理量取得手段6が、物理量記憶部42に記憶された各格子点における物理量φ,ψおよび流束fφ,gψを取得するステップS1。続いて、サブルーチン化された第一領域内解析手段7が、取得された物理量φ,ψおよび流束fφ,gψを用いて所定時間Δt経過後の第一領域内における物理量φik1の数値解を求めるステップS2。続いて、サブルーチン化された第二領域内解析手段8は、取得された物理量φ,ψおよび流束fφ,gψを用いて所定時間Δt後の第二領域内における物理量ψjk1の数値解を求めるステップS3。そして、時間tが所望する時間Tに達したか否かが判別されるステップS4。なお、本第一実施形態では、図3に示すように、第一領域内解析手段7の処理後に第二領域内解析手段8の処理を行っているが、どちらを先に処理してもよく、並列に処理してもよい。また、本第一実施形態では、第一領域内解析手段7および第二領域内解析手段8をサブルーチン化して解析処理を実行しているが、サブルーチン化せずに物理量取得手段6と一連で処理するようにしてもよい。以下、第一領域内解析手段7および第二領域内解析手段8について、より詳細に説明する。図4に示すように、ステップS2における第一領域内解析手段7では、まず、境界流束生成量算出部71が、上記式3により表される流束境界条件式において定義される境界流束生成量γに係る関数γφ,ψに、物理量取得手段6により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界流束生成量γの算出を行うステップS21。算出された境界流束生成量γは、図6に示すように、次のステップの第一領域側境界流束算出部72に送られる。第一領域側境界流束算出部72は、時間方向に離散化された上記式4により表される流束境界条件式に、境界流束生成量算出部71により算出された境界流束生成量γ、および物理量取得手段6により取得された第二領域における境界近傍の流束gψkを代入して、所定時間Δt経過後の第一領域における境界近傍の流束f-φk1を算出するステップS22。算出された第一領域における境界近傍の流束f-φk1は、図6に示すように、次のステップの第一領域内物理量解析部73に送られるとともに、物理量記憶部42に送られる。物理量記憶部42に送られた所定時間Δt経過後の第一領域における境界近傍の流束f-φk1は、第一領域における境界近傍の格子点における流束値として置き換えられる。このように、物理量記憶部42内の物理量を所定時間Δt経過後の物理量に置き換えることにより、時間発展計算を行うことができる。次に、第一領域内物理量解析部73は、時間方向に離散化された上記式5により表される第一領域の輸送方程式に、第一領域側境界流束算出部72により算出された境界近傍の流束f-φk1を代入して導き出される連立方程式を解くことで、第一領域内の各格子点における所定時間Δt経過後の物理量φik1の算出を行うステップS23。算出された第一領域における物理量φik1は、図6に示すように、次のステップの第一領域側境界物理量算出部74に送られるとともに、物理量記憶部42に送られる。第一領域側境界物理量算出部74は、第一領域内物理量解析部73により求められた数値解から所定時間Δt経過後の第一領域における境界近傍の物理量φ-k1を算出するステップS24。そして、算出された第一領域における境界近傍の物理量φ-k1は、図6に示すように、物理量記憶部42に送られる。以上で、第一領域内解析手段7によるステップS2は終了し、次の第二領域内解析手段8によるステップS3へと進む。図5に示すように、ステップS3における第二領域内解析手段8では、まず、境界変化量算出部81が、上記式6の物理量境界条件式により定義された境界変化量cに係る関数cφ,ψに、物理量取得手段6により取得された物理量φおよび物理量ψを代入して、境界における物理量の境界変化量cの算出を行うステップS31。算出された境界変化量cは、図6に示すように、次のステップの第二領域側物理量算出部82に送られる。第二領域側境界物理量算出部82は、時間方向に離散化された上記式7により表される物理量境界条件式に、境界変化量算出部81により算出された境界変化量c、および物理量取得手段6により取得された第一領域における境界近傍の物理量φ-kを代入して、所定時間Δt経過後の第二領域における境界近傍の物理量ψk1を算出するステップS32。算出された第二領域における境界近傍の物理量ψk1は、図6に示すように、次のステップの第二領域内物理量解析部83に送られるとともに、物理量記憶部42に送られる。次に、第二領域内物理量解析部83は、時間方向に離散化された上記式8により表される第二領域の輸送方程式に、第二領域側境界物理量算出部82により算出された境界近傍の物理量ψk1を代入して導き出された連立方程式を解くことによって第二領域内の各格子点における所定時間Δt経過後の物理量ψjk1の算出を行うステップS33。算出された第二領域における物理量ψjk1は、図6に示すように、次のステップの第二領域側境界流束算出部84に送られるとともに、物理量記憶部42に送られる。第一領域側境界流束算出部84は、第二領域内物理量解析部83により求められた数値解から所定時間Δt経過後の第一領域における境界近傍の流束gψk1を算出するステップS34。そして、算出された第二領域における境界近傍の流束gψk1は、図6に示すように、物理量記憶部42に送られる。以上で、第二領域内解析手段8によるステップS3は終了し、次のステップS4へと進む。ステップS4では、時間tが所望する時間Tに達したか否かを行う。時間tが所望する時間Tに達していない場合ステップS4:NO、時間t=tΔtとしてステップS1へもどり、ステップ1からステップS3の処理を繰り返し行う。このように、繰り返し処理を行うことにより第一領域内および第二領域内の各格子点における物理量φ,ψの時間発展計算が行われる。そして、時間tが所望する時間Tに達した場合ステップS4:YES、処理を終了させる。以上のような本第一実施形態の輸送方程式解析プログラム1aおよび輸送方程式解析方法によれば、以下のような効果を得ることができる。1.多層領域間を移動する物理量の数値解析を行うことができる。2.流束境界条件式および物理量境界条件式の2つの境界条件を算出する式を用いることにより、適当な各領域間の境界における物理量を求めることができるため、各領域における輸送方程式の数値解の精度を向上させるとともに、収束性も向上する。3.第一領域側と第二領域側との解析を並列処理することができるため、処理時間を向上させることができる。4.2つの輸送方程式の連成問題であれば物理量の種類については特に限定がなく、応用範囲が広い。次に、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法の第二実施形態について説明する。なお、本第二実施形態のうち、前述した第一実施形態の構成と同一若しくは相当する構成については同一の符号を付して再度の説明を省略する。図7は、本第二実施形態における輸送方程式解析プログラム1bを備えたコンピュータの構成を示すブロック図である。また、図8は、本第二実施形態で解析を行う第一領域、第二領域および第三領域を示す概略図である。図8に示すように、本第二実施形態における第一領域、第二領域および第三領域は、第一領域と第二領域とがこれらの間の第一領域側境界を介して連結されており、第二領域と第三領域とがこれらの間の第三領域側境界を介して連結されている。また、第一領域および第三領域は、それぞれ性質の異なる構造および/または物質からなり、前記各領域の間を移動する所定の物理量は、異なる輸送方程式により表されるものである。具体的には、第一領域における輸送方程式は、上記式1により表されるととともに、第三領域における輸送方程式は、上記式2により表される。さらに、本第二実施形態における第二領域は薄膜状に形成されている。そのため、第二領域内における物理量および流束は、空間的にほぼ一定であると見なすことができる。よって、本第二実施形態では、第二領域内の物理量およびその物理量の流束を、それぞれθおよびhθと近似した。本発明では、以上のような第二領域における物理量等の近似を薄膜近似と定めている。以下、本第二実施形態の各構成について、より詳細に説明する。本第二実施形態におけるコンピュータ1Bは、図7に示すように、主として、入力手段2、表示手段3、記憶手段4および演算処理手段5から構成されている。このうち本第二実施形態における演算処理手段5は、記憶手段4にインストールされた第二実施形態の輸送方程式解析プログラム1bを実行させることにより、図7に示すように、物理量取得手段6、第一領域内解析手段9および第三領域内解析手段10としてコンピュータ1Bを機能させるようになっている。本第二実施形態における第一領域内解析手段9は、上述した第一実施形態における第一領域内解析手段7と同様に、所定時間Δt経過後における第一領域内の物理量の数値解を求めるものである。また、本第二実施形態における第一領域内解析手段9は、図7に示すように、第一領域側境界流束生成量算出部91と、第三領域側境界流束生成量算出部92と、第一領域側境界流束算出部93と、第一領域内物理量解析部94と、第一領域側境界物理量算出部95とを有している。つまり、上述した第一実施形態では、境界が1つしかないため境界の境界流束生成量の算出は、境界流束生成量算出部71による1つの処理で行うことができるのに対し、本第二実施形態では、境界が、第一領域側境界と第三領域側境界の2つあるため、各境界における流束の算出についても第一領域側境界流束生成量算出部91と第三領域側境界流束生成量算出部92とによる2つの処理を要する。第一領域側境界流束生成量算出部91は、以下の式9の流束境界条件式により定義された境界流束生成量γ-に係る関数γ-φ,θを算出するものであり、前記関数γ-φ,θに取得された物理量φおよび物理量θを代入して第一領域側境界における物理量の境界流束生成量γ-を算出している。9また、上記式9で示される流束境界条件式は、第一領域と第二領域との間のノイマン条件による境界条件を表す式であり、第一領域側から第一領域側境界へ流入出する物理量の流束fφと第二領域側から第一領域側境界へ流入出する前記物理量の流束hθとの連続性を示している。ここで、-は、第一領域側境界の値であることを示している。次に、第三領域側境界流束生成量算出部92は、以下の式10の流束境界条件式により定義された境界流束生成量γに係る関数γθ,ψを算出するものであり、前記関数γθ,ψに物理量取得手段6により取得された物理量ψおよびθを代入して第三領域側境界における物理量の境界流束生成量γを算出している。10また、上記式10で示される流束境界条件式は、第三領域と第二領域との間のノイマン条件による境界条件を表す式であり、第三領域側から第三領域側境界へ流入出する物理量の流束gψと第二領域側から第三領域側境界へ流入出する前記物理量の流束hθとの連続性を示している。ここで、は、第三領域側境界の値であることを示している。第一領域側境界流束算出部93は、第一領域側境界流束生成量算出部91により算出された境界流束生成量γ-、第三領域側境界流束生成量算出部92により算出された境界流束生成量γ、および前記物理量取得手段6により取得された第三領域側境界の流束gψkを、前記流束境界条件式に代入して所定時間Δt経過後の第一領域側境界の流束f-φk1を算出するものである。そのため、第一領域側境界流束算出部93では、時間方向に離散化された、以下の式11および式12からなる流束境界条件式が用いられる。1112具体的には、まず、上記式11に第三領域側境界流束生成量算出部92により算出された境界流束生成量γと、第三領域境界の流束gψkを代入して、所定時間Δt経過後の第二領域内の流束hθk1を算出する。そして、上記式12に、第一領域側境界流束生成量算出部91により算出された境界流束生成量γ-と式11により算出された所定時間Δt経過後の第二領域内の流束hθk1を代入して、所定時間Δt経過後の第一領域側境界の流束f-φk1を算出する。本第二実施形態における第一領域内物理量解析部94は、第一実施形態における第一領域内物理量解析部73と同様の機能を有しており、第一領域側境界流束算出部93により算出された第一領域側境界の流束f-φk1を、上記式5で表される第一領域の輸送方程式に代入して所定時間Δt経過後の第一領域内における物理量φik1の数値解を求めるものである。また、本第二実施形態における第一領域側境界物理量算出部95は、第一実施形態における第一領域側境界物理量算出部74と同様の機能を有しており、第一領域内物理量解析部94により求められた数値解から所定時間Δt経過後の第一領域側境界の物理量φ-k1を算出するものである。次に、第三領域内解析手段10について説明する。本第二実施形態における第三領域内解析手段10は、第一実施形態における第二領域内解析手段8と同様に、所定時間Δt経過後における第三領域内の物理量の数値解を求めるものである。本第二実施形態における第二領域内解析手段10は、図7に示すように、第三領域側境界変化量算出部101と、第一領域側境界変化量算出部102と、第三領域側境界物理量算出部103と、第三領域内物理量解析部104と、第三領域側境界流束算出部105とを有している。第三領域側境界変化量算出部101は、以下の式13の流束境界条件式により定義された境界変化量cに係る関数cθ,ψを算出するものであり、前記関数cθ,ψに前記物理量取得手段6により取得された物理量ψおよび物理量θを代入して第三領域側境界における物理量の境界変化量cを算出している。13また、上記式13で示される流束境界条件式は、第三領域と第二領域との間のディレクレ条件による境界条件を表す式であり、第三領域側境界における第三領域側の物理量ψと第二領域側の物理量θとの連続性を示している。第一領域側境界変化量算出部102は、以下の式14の流束境界条件式により定義された境界変化量c-に係る関数c-φ,θを算出するものであり、前記関数c-φ,θに前記物理量取得手段6により取得された物理量φおよび物理量θを代入して第一領域側境界における物理量の境界変化量c-を算出している。14また、上記式14で示される流束境界条件式は、第一領域と第二領域との間のディレクレ条件による境界条件を表す式であり、第一領域側境界における第一領域側の物理量φと第二領域側の物理量θとの連続性を示している。第三領域側境界物理量出部103は、第三領域側境界変化量算出部101により算出された境界変化量c、第一領域側境界変化量算出部102により算出された境界変化量c-、および前記物理量取得手段6により取得された第一領域境界の物理量φ-を、物理量境界条件式に代入して所定時間Δt経過後の第三領域側境界の物理量ψk1を算出するものである。そのため、第三領域側境界物理量算出部103では、時間方向に離散化された、以下の式15および式16からなる物理量境界条件式が用いられる。1516具体的には、まず、上記式15に第一領域側境界変化量算出部102により算出された境界変化量c-と、第一領域境界の物理量φ-を代入して、所定時間Δt経過後の第二領域内の物理量θk1を算出する。そして、上記式16に、第三領域側境界変化量算出部101により算出された境界変化量cと、式15により算出された所定時間Δt経過後の第二領域内の物理量θk1を代入して、所定時間Δt経過後の第三領域側境界の物理量ψ-k1を算出する。本第二実施形態における第三領域内物理量解析部104は、第一実施形態における第二領域内物理量解析部83と同様の機能を有しており、第三領域側境界物理量算出部103により算出された第三領域側境界の物理量ψ-k1を、上記式8で表される第三領域の輸送方程式に代入して所定時間Δt経過後の第一領域内における物理量φik1の数値解を求めるものである。また、本第二実施形態における第三領域側境界流束算出部105は、第一実施形態における第二領域側境界流束算出部84と同様の機能を有しており、第三領域内物理量解析部104により求められた数値解から所定時間Δt経過後の第三領域側境界の物理量の流束gψk1を算出するものである。次に、本第二実施形態の輸送方程式解析プログラム1bにおける各構成の作用および輸送方程式解析方法について、図9ないし図11に示すフローチャート図および図12に示す各構成の機能を示す機能ブロック図を用いて説明する。図9および図12に示すように、本第二実施形態の輸送方程式解析プログラム1bおよび輸送方程式解析方法では、まず、物理量取得手段6が、物理量記憶部42に記憶された各格子点における物理量φ,ψおよび流束fφ,gψを取得するステップSS1。続いて、サブルーチン化された第一領域内解析手段9が、取得された物理量φ,ψおよび流束fφ,gψを用いて所定時間Δt経過後の第一領域内における物理量φik1の数値解を求めるステップSS2。続いて、サブルーチン化された第三領域内解析手段10は、取得された物理量φ,ψおよび流束fφ,gψを用いて所定時間Δt経過後の第三領域内における物理量ψjk1の数値解を求めるステップSS3。そして、時間tが所望する時間Tに達したか否かが判別されるステップSS4。なお、第一領域内解析手段9の処理と第三領域内解析手段10の処理とはどちらを先に処理してもよく、並列に処理してもよい。また、本第二実施形態では、第一領域内解析手段9および第三領域内解析手段10をサブルーチンとせずに物理量取得手段6と一連で処理するようにしてもよい。図10に示すように、ステップSS2における第一領域内解析手段9では、まず、第一領域側境界流束生成量算出部91が、上記式9の第一領域側境界における流束境界条件式により定義された境界流束生成量γ-に係る関数γ-φ,θに、物理量取得手段6により取得された物理量φおよび物理量θを代入して、第一領域側境界における物理量の境界流束生成量γ-の算出を行うステップSS21。算出された境界流束生成量γ-は、図12に示すように、第一領域側境界流束算出部93に送られる。また、第三領域側境界流束生成量算出部92が、上記式10により表される第三領域側境界における流束境界条件式により定義され境界流束生成量γに係る関数γθ,ψに、物理量取得手段6で取得された物理量ψおよびθを代入して、第一領域側境界における物理量の境界流束生成量γの算出を行うステップSS22。算出された境界流束生成量γは、図12に示すように、第一領域側境界流束算出部93に送られる。第一領域側境界流束算出部93は、時間方向に離散化された上記式11により表される第三領域側境界における流束境界条件式に、第三領域側境界流束生成量算出部92により算出された境界流束生成量γ、および物理量取得手段6により取得された第三領域側境界の流束gψkを代入して、所定時間Δt経過後の第二領域内の流束hθk1を算出するとともに、上記式12に、第一領域側境界流束生成量算出部91により算出された境界流束生成量γ-と式11により算出された所定時間Δt経過後の第二領域内の流束hθk1を代入して、所定時間Δt経過後の第一領域側境界の流束f-φk1を算出するステップSS23。算出された第一領域側境界の流束f-φk1は、図12に示すように、次のステップの第一領域内物理量解析部94に送られるとともに、物理量記憶部42に送られる。次に、第一領域内物理量解析部94は、時間方向に離散化された上記式5により表される第一領域の輸送方程式に、第一領域側境界流束算出部93により算出された第一領域側境界の流束f-φk1を代入して導き出される連立方程式を解くことにより、第一領域内の各格子点における所定時間Δt経過後の物理量φik1の算出を行うステップSS24。算出された第一領域における物理量φik1は、図12に示すように、次のステップの第一領域側境界物理量算出部95に送られるとともに、物理量記憶部42に送られる。第一領域側境界物理量算出部95は、第一領域内物理量解析部94により求められた数値解から所定時間Δt経過後の第一領域における境界近傍の物理量φ-k1を算出するステップSS25。そして、算出された第一領域における境界近傍の物理量φ-k1は、図12に示すように、物理量記憶部42に送られる。図11に示すように、ステップSS3における第二領域内解析手段10では、まず、第三領域側境界変化量算出部101が、上記式13の第三領域側境界における物理量境界条件式により定義された境界変化量cθ,ψに係る関数に、物理量取得手段6により取得された物理量ψおよび物理量θを代入して、第三領域側境界における物理量の境界変化量cの算出を行うステップSS31。算出された境界変化量cは、図11に示すように、第三領域側物理量算出部103に送られる。次に、第一領域側境界変化量算出部102が、上記式14の第一領域側境界における物理量境界条件式により定義された境界変化量c-φ,θに係る関数に、物理量取得手段6により取得された物理量φおよび物理量θを代入して、第一領域側境界における物理量の境界境界変化量c-の算出を行うステップSS32。算出された境界変化量c-は、図11に示すように、第三領域側物理量算出部103に送られる。第三領域側境界物理量算出部103は、時間方向に離散化された上記式15により表される第一領域側境界における流束境界条件式に、第一領域側境界変化量算出部102により算出された境界変化量c-、および物理量取得手段6により取得された第一領域側境界の物理量φkを代入して、所定時間Δt経過後の第二領域内の物理量θk1を算出するとともに、上記式16に、第三領域側境界変化量算出部101により算出された境界変化量cと、式15により算出された所定時間Δt経過後の第二領域内の物理量θk1を代入して、所定時間Δt経過後の第三領域側境界の物理量ψk1を算出するステップSS23。算出された第三領域側境界の物理量ψk1は、図12に示すように、次のステップの第三領域内物理量解析部104に送られるとともに、物理量記憶部42に送られる。次に、第三領域内物理量解析部104は、時間方向に離散化された上記式8により表される第三領域の輸送方程式に、第三領域側境界物理量算出部103により算出された第三領域側境界における物理量ψk1を代入して導き出される連立方程式を解くことにより、第三領域内の各格子点における所定時間Δt経過後の物理量ψjk1の算出を行うステップSS34。算出された第三領域側境界における物理量ψjk1は、図12に示すように、次のステップの第三領域側境界流束算出部105に送られるとともに、物理量記憶部42に送られる。第三領域側境界流束算出部105は、第三領域内物理量解析部104により求められた数値解から所定時間Δt経過後の第三領域側境界における流束gψk1を算出するステップSS35。そして、算出された第三領域側境界の流束gψk1は、図12に示すように、物理量記憶部42に送られる。ステップSS4では、時間tが所望する時間Tに達したか否かを行う。時間tが所望する時間Tに達していない場合ステッSプS4:NO、時間t=tΔtとし、ステップSS1からステップSS3の処理を繰り返し行う。そして、時間tが所望する時間Tに達した場合ステップSS4:YES、処理を終了させる。以上のような本第二実施形態の輸送方程式解析プログラム1bおよび輸送方程式解析方法によれば、薄膜近似した第二領域にける格子解像度を高くする必要が無くなるため、計算負荷を小さくすることができるという効果が得られる。本発明を用いた高分子形燃料電池の解析本実施例1では、上述した第一実施形態の解析方法を用いて、高分子形燃料電池内における各種物理量の輸送方程式および各層の境界における境界条件式を例示しつつ、解析について説明する。高分子形燃料電池は、図13に示すように、電解膜と、この電界膜を中心に陰極側および陽極側にそれぞれ触媒層と多孔質層とが多層状に連結されている。また、高分子形燃料電池では、水素や酸素などの作動ガス、生成物である水または水蒸気等が混合された状態で前記各層内を移動する、いわゆる複合流動場が構成されている。そこで、本実施例1では、各層における物理量を液水s、水蒸気質量分率Ywおよび電界質内の含水量λに分けて、それぞれについて輸送方程式が設定されている。また、液水sおよび水蒸気質量分率Ywは直接相変化が行われるとともに、水蒸気質量分率Ywおよび電界質内の含水量λは直接相変化が行われるものとし、一方、液水sおよび電界質内の含水量λは直接相変化が行われないものとした。1各層における輸送方程式について各層における輸送方程式について、図14を参照しながら説明する。図14に示すように、以下の式17は、触媒層内における液水sの輸送方程式である。17また、以下の式18は、多孔質層内における液水sの輸送方程式である。18ここで、上記式17および上記式18内の各係数は、以下の式19ないし式21の関係を有している。192021また、εは多効率、ρlは水の密度、ρgはガスの密度、u気体の速度、Dc拡散係数、κ透過率、νは動粘性計数、gは重力加速度、Mlはモル質量、kcは凝縮係数、keは蒸発係数、Xwはモル分率、Rは気体定数、T温度、P圧力、Psatは飽和水蒸気量、σは表面張力、θcは液水sの接触角およびμは粘性係数を示している。また、m1は液水sと水蒸気質量分率Ywの間の相変化量である。なお、電界膜は固体であり液水sは入り込めないため、図14に示すように、電界膜内の液水s=0である。次に、水蒸気質量分率Ywの輸送方程式について説明する。以下の式22は、触媒層内における水蒸気質量分率Ywの輸送方程式である。22また、以下の式23は、多孔質層内における水蒸気質量分率Ywの輸送方程式である。23また、上記式22および上記式23内の各係数は、以下の式24の関係を有している。24ここで、Yiは化学種iの質量分率、Di,effは参照拡散係数、Di0は拡散係数、P0は基準圧力、T0は基準温度、ωは化学反応による水または水蒸気の発生量、m2は水蒸気質量分率Ywと電界質内の含水量λと間の相変化量である。また、水蒸気質量分率Ywは、液水sと同様に、電界膜内に入り込めないため、電界膜内では水蒸気質量分率Yw=0である。次に、電界質内の含水量λの輸送方程式について説明する。以下の式25は、電解膜における電界質内の含水量λの輸送方程式である。25また、以下の式26は、触媒層内における電解質内の含水量λの輸送方程式である。26また、上記式25および上記式26内の各係数は、以下の式27および式28の関係を有している。2728ここで、ρmは電解膜の密度、EWは電解膜の等価質量、Vexは電解膜の膨張係数、Dwlは電解膜の拡散係数、Diは拡散係数、ndは電気浸透抵抗係数、iは電流、Fはファラデー定数である。なお、電界膜は固体であり液水sは入り込めないため、電界膜内では液水s=0である。以上のように、高分子形燃料電池の問題においては、液水sは触媒層と多孔質層の二つの領域を移動する物理量であり、各領域にいて異なる輸送方程式により示すことができる。また、水蒸気質量分率Ywは触媒層と多孔質層の二つの領域を移動する物理量であり、各領域において異なる輸送方程式により示すことができる。さらに、電界質内の含水量λは電解質と触媒層の二つの領域を移動する物理量であり、各領域において異なる輸送方程式により示すことができる。すなわち、高分子形燃料電池内の液水s、水蒸気質量分率Ywおよび電界質内の含水量λは、それぞれ性質の異なる第一領域および第二領域が第一領域と第二領域との間の境界を介して連結されており、各々の領域を移動する物理量が第一領域および第二領域で異なる輸送方程式によって表されるものであるため、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法によって処理可能な問題である。2液水sの輸送方程式の解析についてまず、液水sの輸送方程式の解析について説明する。本実施例1では、第一領域を触媒層、第二領域を多孔質層とした。液水sの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式29で表される。29また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式30で表される。30ただし、第一領域側のPc,CLは、以下の式31および式32によって表されるものである。3132なお、Pcはsの三次の多項式であり、逆関数Pc-1を直接求めることが難しい。そこで、以下の方法でPcの線形化を行った。まず、s=aにおいてテイラー展開を行うと、以下の式33で表される一次多項式が得られる。33上記式33をPcに適用すると、以下の式34が得られる。34以上を踏まえると、上記式31は線形化され、以下の式35として表される。35また、第二領域側のPc,GDLにも適用すると、Pc,GDLは以下の式36として表される。36ただし、上記式35および上記式36は以下の関係を有する。ここで、aは前時刻のsa=sk-1を用いることとすると、上記式30は、以下の式37として表される。37第一領域側の輸送方程式である上記式17を離散化すると、以下の式38として表される。38また、第二領域側の輸送方程式である上記式18を離散化すると、以下の式39として表される。39なお、上記式38および上記式39は、時間発展計算をオイラー陽解法で解く場合の式である。次に、各境界条件についても離散化を行う。以下の式40は、ノイマン条件における上記式29を離散化したものである。40また、以下の式41は、ディレクレ条件における上記式37を離散化したものである。41次に、以上の各式を用いて、第一領域および第二領域における液水sの輸送方程式の解析について、特に、第一領域内解析手段7および第二領域内解析手段8の解析手順に従って、以下に説明する。第一領域内解析手段7では、まず、境界流束生成量算出部71が、上記式29により境界流束生成量γを算出する。本件の場合は、図14に示すように、境界における流束Jfluxは、化学反応等により変化しないと仮定しているため、境界流束生成量γ=0である。次に、第一領域側境界流束算出部72は、上記式40を用いて第一領域の境界近傍の流束f-sCLk1を算出する。なお、本件における流束は、各領域における多孔率εに影響されるため、α=εCLであり、β=εGDLである。第一領域内物理量解析部73は、第一領域側境界流束算出部72より算出された流束f-sCLk1を用いて上記式38を解き、液水sCLk1を求める。そして、第一領域側境界物理量算出部74は、第一領域内物理量解析部73により求められた液水sCLk1の値を用いて、第一領域の境界近傍における液水sの値s-k1を算出する。以上により、第一領域内の所定時間Δt経過後における液水sの値を全て求めることができる。次に、第二領域内解析手段8では、まず、境界変化量算出部81が、上記式37により境界における液水sの境界変化量cを算出する。本件では、境界におけるPcの値は化学反応等により変化しないとしているため、境界変化量c=0である。第二領域側境界物質量算出部82は、上記式41を用いて第二領域の境界近傍の液水sCLk1を算出する。なお、本件では、a=1でありb=1である。第二領域内物理量解析部83は、第二領域側境界物理量算出部82より算出された液水sCLk1を用いて上記式39を解き、液水GDLk1を求める。そして、第二領域側境界流束算出部は、第二領域内物理量解析部83により求められた液水sGDLk1の値を用いて、第二領域の境界近傍における液水sの流束gsGDLk1を算出する。以上により、第二領域内の所定時間Δt経過後における液水sの値を全て求めることができる。3水蒸気質量分率Ywの輸送方程式の解析について次に、水蒸気質量分率Ywの輸送方程式の解析について説明する。本実施例1では、第一領域を触媒層、第二領域を多孔質層とした。水蒸気質量分率Ywの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式42で表される。42また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式43で表される。43次に、第一領域側の輸送方程式である上記式22を離散化すると、以下の式44として表される。44また、第二領域側の輸送方程式である上記式23を離散化すると、以下の式45として表される。45なお、上記式44および上記式45は、時間発展計算をオイラー陽解法で解く場合の式である。次に、各境界条件についても離散化を行う。以下の式46は、ノイマン条件における上記式42を離散化したものである。46また、以下の式47は、ディレクレ条件における上記式43を離散化したものである。47次に、以上の各式を用いて、第一領域および第二領域における水蒸気質量分率Ywの輸送方程式の解析について、特に、第一領域内解析手段7および第二領域内解析手段8の解析手順に従って、以下に説明する。第一領域内解析手段7では、まず、境界流束生成量算出部71が、上記式42により境界流束生成量γを算出する。本件の場合は、図14に示すように、境界における流束Jfluxは、化学反応等により変化しないと仮定しているため、境界流束生成量γ=0である。次に、第一領域側境界流束算出部72は、上記式46を用いて第一領域の境界近傍の流束f-YCLk1を算出する。なお、本件における流束は、各領域における多孔率εに影響されるため、α=εCLであり、β=εGDLである。第一領域内物理量解析部73は、第一領域側境界流束算出部72より算出された流束f-sCLk1を用いて上記式44を解き、水蒸気質量分率YCLk1を求める。そして、第一領域側境界物理量算出部74は、第一領域内物理量解析部73により求められた水蒸気質量分率YCLk1の値を用いて、第一領域の境界近傍における水蒸気質量分率Yの値Y-k1を算出する。以上により、第一領域内の所定時間Δt経過後における水蒸気質量分率Ywの値を全て求めることができる。次に、第二領域内解析手段8では、まず、境界変化量算出部81が、上記式43により境界における水蒸気質量分率Ywの境界変化量cを算出する。本件では、境界におけるYwの値は化学反応等により変化しないと仮定しているため、境界変化量c=0である。第二領域側境界物質量算出部82は、上記式47を用いて第二領域の境界近傍の水蒸気質量分率YCLk1を算出する。なお、本件では、a=1でありb=1である。第二領域内物理量解析部83は、第二領域側境界物理量算出部82より算出された水蒸気質量分率YCLk1を用いて上記式45を解き、水蒸気質量分率YGDLk1を求める。そして、第二領域側境界流束算出部は、第二領域内物理量解析部83により求められた水蒸気質量分率YGDLk1の値を用いて、第二領域の境界近傍における水蒸気質量分率Ywの流束gYGDLk1を算出する。以上により、第二領域内の所定時間Δt経過後における水蒸気質量分率Ywの値を全て求めることができる。4電界質内の含水量λの輸送方程式の解析について次に、電界質内の含水量λの輸送方程式の解析について説明する。本実施例1では、第一領域を電解膜、第二領域を触媒層とした。電界質内の含水量λの輸送方程式を解析するための、第一領域と第二領域との間の境界におけるノイマン条件は、以下の式48で表される。48また、第一領域と第二領域との間の境界におけるディレクレ条件は、以下の式49で表される。49次に、第一領域側の輸送方程式である上記式25を離散化すると、以下の式50として表される。50また、第二領域側の輸送方程式である上記式26を離散化すると、以下の式51として表される。51なお、上記式50および上記式51は、時間発展計算をオイラー陽解法で解く場合の式である。次に、各境界条件についても離散化を行う。以下の式52は、ノイマン条件における上記式48を離散化したものである。52また、以下の式53は、ディレクレ条件における上記式49を離散化したものである。53次に、以上の各式を用いて、第一領域および第二領域における電界質内の含水量λの輸送方程式の解析について、特に、第一領域内解析手段7および第二領域内解析手段8の解析手順に従って、以下に説明する。第一領域内解析手段7では、まず、境界流束生成量算出部71が、上記式48により境界流束生成量γを算出する。本件の場合は、図14に示すように、境界における流束Jfluxは、化学反応等により変化しないと仮定しているため、境界流束生成量γ=0である。次に、第一領域側境界流束算出部72は、上記式46を用いて第一領域の境界近傍の流束f-λCLk1を算出する。なお、本件では、α=1であり、β=1である。第一領域内物理量解析部73は、第一領域側境界流束算出部72より算出された流束f-λCLk1を用いて上記式52を解き、電界質内の含水量λCLk1を求める。そして、第一領域側境界物理量算出部74は、第一領域内物理量解析部73により求められた電界質内の含水量λCLk1の値を用いて、第一領域の境界近傍における電界質内の含水量λの値λ-k1を算出する。以上により、第一領域内の所定時間Δt経過後における電界質内の含水量λの値を全て求めることができる。次に、第二領域内解析手段8では、まず、境界変化量算出部81が、上記式49により境界における電界質内の含水量λの境界変化量cを算出する。本件では、境界におけるλの値は化学反応等により変化しないと仮定しているため、境界変化量c=0である。第二領域側境界物質量算出部82は、上記式53を用いて第二領域の境界近傍の電界質内の含水量λCLk1を算出する。なお、本件では、a=1でありb=1である。第二領域内物理量解析部83は、第二領域側境界物理量算出部82より算出された電界質内の含水量λCLk1を用いて上記式51を解き、電界質内の含水量λGDLk1を求める。そして、第二領域側境界流束算出部84は、第二領域内物理量解析部83により求められた電界質内の含水量λGDLk1の値を用いて、第二領域の境界近傍における電界質内の含水量λの流束gλGDLk1を算出する。以上により、第二領域内の所定時間Δt経過後における電界質内の含水量λwの値を全て求めることができる。本発明における薄膜近似を用いた高分子形燃料電池の解析本実施例2では、上述した第二実施形態に基づく、本発明における薄膜近似を用いた高分子形燃料電池の解析手法について説明する。本実施例2における解析場は、図15に示すように、電解膜、触媒層および多孔質層の3つの領域からなり、中間に位置する触媒層における物理量を薄膜近似して処理するものである。1各層における輸送方程式について各層における輸送方程式について、図15を参照しながら説明する。図15に示すように、電解膜および多孔質層における液水s、水蒸気質量分率Ywおよび電解質内の含水量λの輸送方程式は、上述した実施例1におけるものと同一である。具体的には、電解膜内では、液水s=0、水蒸気質量分率Yw=0であるとともに、電解質内の含水量λの輸送方程式は上記式25により表される。また、多孔質層内では、液水sの輸送方程式は上記式18により表され、水蒸気質量分率Ywの輸送方程式は上記式23により表されるとともに、電解質内の含水量λ=0である。なお、以下の説明においては、上記式25は、基礎方程式形で表された、以下の式54として表した。54同様に、上記式18は、以下の式55とした。55さらに、上記式23は、以下の式56とした。56また、本実施例2では、触媒層内の物理量の輸送は一定であると仮定して、その平均値s、Yw、λおよび平均流束hs、hYw、hλとして定義した。2各輸送方程式の解析について本実施例2では、液水s、水蒸気質量分率Ywおよび電解質内の含水量λを同時に処理する。また、多孔質層を第一領域、触媒層を第二領域および電解膜を第三領域とした。電解膜から多孔質層への物理量の移動に対しては、ノイマン条件を与えた。ただし、電解膜における電解質内の含水量λは、実際には多孔質層内の水蒸気質量分率Ywおよび液水sとして移動するが、すべてが多孔質層内の水蒸気質量分率Ywとして移動するものとして仮定した。なお、この仮定下においては、多孔質層内の水蒸気質量分率Ywが飽和量を超えた際には直ちに液水sに変換される関係にあり、電解質内の含水量λは間接的に液水sになるものであるため、物理的な関連性は保たれているものと考えられる。以上の仮定を踏まえると、第三領域側境界におけるノイマン条件は、以下の式57で表される。57また、第一領域側境界におけるノイマン条件は、以下の式58で表される。58また、多孔質層から電解膜への物理量の移動に対しては、ディレクレ条件を与えた。ここで、触媒層内における電解質内の含水量λは、平衡状態に達していると考え、以下の式59に示すように、液水sとawwateractivityの関数で表されるものとした。59第三領域側境界におけるディレクレ条件は、以下の式60で表される。60また、第一領域境界においては、液水sおよび水蒸気質量分率Ywのディレクレ条件が与えられ、液水sに関しては、以下の式61で表される。61また、水蒸気質量分率Ywに関しては、以下の式62で表される。62次に、上述した各式の離散化を行う。第一領域における電解質内の含水量λの輸送方程式である上記式54の離散化をすると、以下の式63として表される。63また、第三領域における液水sの輸送方程式である上記式55の離散化をすると、以下の式64として表される。64さらに、第三領域における水蒸気質量分率Ywの輸送方程式である上記式56を離散化すると、以下の式65として表される。65次に、各境界条件についても離散化を行う。第三領域側境界におけるノイマン条件の上記式57を離散化すると、以下の式66として表される。66また、本実施例2では、電解質内の含水量λは、すべてが多孔質層内の水蒸気質量分率Ywとして移動するものと仮定したため、電解質内の含水量λと水蒸気質量分率Ywとの間には、以下の式67の関係が成り立つ。67また、第一領域側境界におけるノイマン条件の上記式58の離散化を行うと、以下の式68として表される。68次に、第三領域側境界のディレクレ条件の上記式60を離散化すると、以下の式69として表される。69また、第三領域側境界における液水sのディレクレ条件の上記式61の離散化をすると、以下の式70として表される。70さらに、第一領域側境界における水蒸気質量分率Ywのディレクレ条件の上記式62を離散化すると、以下の式71として表される。71また、電解質内の含水量λを表した上記式59を離散化すると、以下の式72として表される。72次に、以上の各式を用いて、第一領域および第三領域における各物理量の輸送方程式の解析について、特に、第二実施形態における第一領域内解析手段9および第三領域内解析手段10の解析手順に従って、以下に説明する。第一領域内解析手段9では、まず、第一領域側境界流束生成量算出部91が、境界流束生成量γを算出する。本件では、境界流束生成量γは、以下の式73として算出される。73なお、陰極側の計算には、以下の式74で表されるバトラーボルマー方程式を用いた。74ここで、電子2mol流れたときにH2Oは1mol発生するため、以下の式75で表される水が発生する。75この式75に、触媒層の厚みLを掛けると、上記式73の陰極側の値が導き出される。次に、第三領域側境界境界流束生成量算出部92が、境界流束生成量γ-を算出する。本件では、境界流束生成量γ-=0である。次に、第一領域側境界流束算出部93は、上記式66、上記式67および上記式68を用いて第一領域側境界の流束gYw、GDLk1を算出する。なお、本件におけるαおよびβを以下に示す。第一領域内物理量解析部94は、第一領域側境界流束算出部73より算出された流束gYw、GDLk1を用いて上記式65を解き、多孔質層内の水蒸気質量分率Yjk1を求める。そして、第一領域側境界物理量算出部95は、第一領域内物理量解析部94により求められた多孔質層内の水蒸気質量分率Yjk1の値を用いて、第一領域側境界における水蒸気質量分率Ywの値Yk1を算出する。以上により、第一領域内の所定時間Δt経過後における水蒸気質量分率Ywの値を全て求めることができる。次に、第三領域内解析手段10では、第三領域側境界変化量算出部101が、上記式60により境界変化量c-を算出する。本件では、境界変化量c-=0である。また、第一領域側境界変化量算出部102が、上記式61および上記式62により、境界変化量cs,および境界変化量cY,を算出する。本件では、境界変化量cs,=0であり、境界変化量cY,=0である。次に、第三領域側境界物質量算出部103は、上記式69、上記式71および上記式72を用いて第三領域側境界の電界質内の含水量λ-k1を算出する。なお、本件では、aおよびbに関する値はすべて1である。第三領域内物理量解析部104は、第三領域側境界物質量算出部103より算出された電界質内の含水量λ-k1を用いて上記式63を解き、電界質内の含水量λik1を求める。そして、第三領域側境界流束算出部105は、第三領域内物理量解析部104により求められた電界質内の含水量λik1の値を用いて、第三領域側境界における電界質内の含水量λの流束fλGDLk1を算出する。以上により、第三領域内の所定時間Δt経過後における電界質内の含水量λwの値を全て求めることができる。なお、本発明に係る輸送方程式解析プログラムおよび輸送方程式解析方法による解析対象は、前述した高分子燃料電池の問題に限定されるものではなく、他にも適用することができる。例えば、二次電池におけるイオン伝導問題や多孔質内における燃焼問題等に使用することができる。"