(57)【要約】【課題】気液二相流れ評価シミュレーションにおいて、多孔体の気孔率分布が複数の気孔率ピークを有するときに、実験結果との適合性のよい計算を行うことである。【解決手段】多孔体の気孔率分布が複数の気孔率ピークを有するとき、多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成されるものとしてモデル化し（Ｓ１０）、各層ごとに、その気孔率分布に応じて、気体相対透過率ｋGと液体相対透過率ｋLをそれぞれ算出し（Ｓ１２）、第１層に流入圧力ＰINを与えて、第１層の流出圧力側のＰOUT1を求め（Ｓ１４）、第１層の流出圧力ＰOUT1を第２層に対する流入圧力ＰIN2として与え、第２層の流出圧力ＰOUT2を求め（Ｓ１６）、これを繰り返して、流出側層について求められた流出圧力をＰOUTとして多孔体の圧力損失ΔＰを算出する（Ｓ２０）。【選択図】図３
2008276357.tif 000002
【特許請求の範囲】
【請求項１】
多孔体における気液二相流れについて、多孔体の気孔率と、多孔体の構造で定まる絶対透過率Ｋと、多孔体における気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ与えて、所定の関係式に基づいて、多孔体の流入側と流出側との間の差圧である圧力損失ΔＰを求める気液二相流れ評価シミュレーション装置であって、
多孔体の気孔率分布が複数の気孔率ピークを有するとき、
多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとしてモデル化する手段と、
各層ごとに、その気孔率分布に応じて、絶対透過率Ｋと、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ割り当てる手段と、
積層モデルの各層のうち、気液二相流れが多孔体に供給される流入側層を第１層とし、流入圧力ＰINを与えて、第１層の圧力損失ΔＰ1から第１層の流出圧力側のＰOUT1を求める手段と、
積層モデルの各層のうち、第１層の次に流出側に配置された層を第２層とし、第１層の流出圧力ＰOUT1を第２層に対する流入圧力ＰIN2として与え、第２層の圧力損失ΔＰ2から第２層の流出圧力ＰOUT2を求める手段と、
積層モデルを構成する層数に応じてこれを繰り返し、積層モデルの各層のうち、気液二相流れが多孔体から流出する流出側層について求められた流出圧力をＰOUTとし、
ΔＰ＝ＰIN−ＰOUT
を多孔体の圧力損失ΔＰとして出力する手段と、
を備えることを特徴とする気液二相流れ評価シミュレーション装置。
【請求項２】
気体密度ρGと気体粘性係数μGと気体質量速度ＧGを有する気体と、液体密度ρLと気体粘性係数μGと液体質量速度ＧLを有する液体とが気液二相流れとして、厚みがＬである多孔体に供給されるとき、
入力条件として、多孔体の気孔率と、多孔体の構造で定まる絶対透過率Ｋと、を与え、
気孔率に基づいて液体の飽和度ｓと毛管圧ｐCを求め、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式に基づいて、気体と液体のそれぞれについて、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ求め、
ΔＰ／Ｌ＝（μGＧG）／（ｋGＫρG）＝（μLＧL）／（ｋLＫρL）
の関係式に基づき圧力損失ΔＰを求める気液二相流れ評価シミュレーション装置であって、
多孔体の気孔率分布が複数の気孔率ピークを有するとき、
多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとし、各層をｉ層として、ｉ層のそれぞれについて、その気孔率分布に応じて、厚さＬiと絶対透過率Ｋiとを割り当てるモデル化手段と、
ｉ層のそれぞれについて、その気孔率分布に応じて、気体相対透過率ｋGiと液体相対透過率ｋLiとを求める相対透過率算出手段と、
積層モデルの各層のうち、気液二相流れが多孔体に供給される流入側層を第１層とし、流入圧力ＰINを与え、ｉ＝１として、
ΔＰ／Ｌ＝（ＰIN−ＰOUT1）／Ｌ1＝（μGＧG）／（ｋG1Ｋ1ρG）＝
（μLＧL）／（ｋL1Ｋ1ρL）
から、第１層の流出圧力ＰOUT1を求める手段と、
積層モデルの各層のうち、第１層の次に流出側に配置された層を第２層とし、第１層の流出圧力ＰOUT1を第２層に対する流入圧力ＰIN2として与える手段と、
第２層について、流入圧力ＰIN2＝ＰOUT1を与え、ｉ＝２として、
ΔＰ／Ｌ＝（ＰIN2−ＰOUT2）／Ｌ2＝（μGＧG）／（ｋG2Ｋ2ρG）＝
（μLＧL）／（ｋL2Ｋ2ρL
から、第２層の流出圧力ＰOUT2を求める手段と、
積層モデルを構成する層数に応じてこれを繰り返し、積層モデルの各層のうち、気液二相流れが多孔体から流出する流出側層について求められた流出圧力をＰOUTとし、
ΔＰ＝ＰIN−ＰOUT
を多孔体の圧力損失ΔＰとして出力する手段と、
を備えることを特徴とする気液二相流れ評価シミュレーション装置。
【請求項３】
請求項２に記載の気液二相流れ評価シミュレーション装置において、
相対透過率算出手段は、
各ｉ層について、ｉ層の気孔率分布を構成する各気孔径ＤPについて液体の飽和度ｓiを割り当てて飽和度分布ｓi＝ｆi（Ｄp）を求め、この飽和度分布ｓiに対応する毛管圧分布ｐCi＝ｇi（ｓ）を求める手段と、
各ｉ層についての飽和度分布ｓi＝ｆi（Ｄp）と毛管圧分布ｐCi＝ｇi（ｓ）とに基づき、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式から、各ｉ層についての気体相対透過率の分布ｋG＝ｋG（ｓ）と液体相対透過率の分布ｋL＝ｋL（ｓ）を求める手段と、
を含むことを特徴とする気液二相流れ評価シミュレーション装置。
【請求項４】
請求項３に記載の気液二相流れ評価シミュレーション装置において、
ｉ層を複数の分割要素に分割し、ｉ層の気体相対透過率の分布ｋG＝ｋG（ｓ）と液体相対透過率の分布ｋL＝ｋL（ｓ）に基づいて、各分割要素を分割要素Ｊとして、分割要素Ｊのそれぞれに気体相対透過率ｋGJと液体相対透過率ｋLJを割り当て、有限要素法を用いて多孔体の圧力損失ΔＰを求めることを特徴とする気液二相流れ評価シミュレーション装置。
【請求項５】
気体密度ρGと気体粘性係数μGと気体質量速度ＧGを有する気体と、液体密度ρLと気体粘性係数μGと液体質量速度ＧLを有する液体とが気液二相流れとして、厚みがＬである多孔体に供給されるとき、
入力条件として、多孔体の気孔率と、多孔体の構造で定まる絶対透過率Ｋと、を与え、
気孔率に基づいて液体の飽和度ｓと毛管圧ｐCを求め、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式に基づいて、気体と液体のそれぞれについて、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ求め、
ΔＰ／Ｌ＝（μGＧG）／（ｋGＫρG）＝（μLＧL）／（ｋLＫρL）
の関係式に基づき圧力損失ΔＰを求める気液二相流れ評価シミュレーション装置において実行される気液二相流れ評価シミュレーションプログラムであって、
多孔体の気孔率分布が複数の気孔率ピークを有するとき、
多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとし、各層をｉ層として、ｉ層のそれぞれについて、その気孔率分布に応じて、厚さＬiと絶対透過率Ｋiとを割り当てるモデル化処理手順と、
ｉ層のそれぞれについて、その気孔率分布に応じて、気体相対透過率ｋGiと液体相対透過率ｋLiとを求める相対透過率算出処理手順と、
積層モデルの各層のうち、気液二相流れが多孔体に供給される流入側層を第１層とし、流入圧力ＰINを与え、ｉ＝１として、
ΔＰ／Ｌ＝（ＰIN−ＰOUT1）／Ｌ1＝（μGＧG）／（ｋG1Ｋ1ρG）＝
（μLＧL）／（ｋL1Ｋ1ρL）
から、第１層の流出圧力ＰOUT1を求める処理手順と、
積層モデルの各層のうち、第１層の次に流出側に配置された層を第２層とし、第１層の流出圧力ＰOUT1を第２層に対する流入圧力ＰIN2として与える処理手順と、
第２層について、流入圧力ＰIN2＝ＰOUT1を与え、ｉ＝２として、
ΔＰ／Ｌ＝（ＰIN2−ＰOUT2）／Ｌ2＝（μGＧG）／（ｋG2Ｋ2ρG）＝
（μLＧL）／（ｋL2Ｋ2ρL）
から、第２層の流出圧力ＰOUT2を求める処理手順と、
積層モデルを構成する層数に応じてこれを繰り返し、積層モデルの各層のうち、気液二相流れが多孔体から流出する流出側層について求められた流出圧力をＰOUTとし、ΔＰ＝ＰIN−ＰOUTを多孔体の圧力損失ΔＰとして出力する処理手順と、
が実行される気液二相流れ評価シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、気液二相流れ評価シミュレーション装置及び気液二相流れ評価シミュレーションプログラムに係り、特に、多孔体の気孔率分布が複数の気孔率ピークを有するときの気液二相流れ評価シミュレーション装置及び気液二相流れ評価シミュレーションプログラムに関する。
【背景技術】
【０００２】
気孔を有するポーラスな物質である多孔体は、気体や液体を透過するので、ろ過層や拡散層等として用いられる。このようにろ過層や拡散層等として用いられるものについては、気体や液体がどのように透過するかを評価する必要がある。例えば、多孔体の一方側を供給側として気体や液体を供給し、他方側を排出側として透過した気体や液体を排出させ、そのときの供給側の圧力と排出側の圧力との差圧を、多孔体の圧力損失として評価することが行われる。
【０００３】
多孔体は固体であって気孔を有するものであるので、気体や液体がどのように透過するかは、気孔率や気孔における毛管圧、含水率等の多孔体の特性に左右される。したがって、多孔体の圧力損失を予めシミュレーション等で求めようとするときには、これら多孔体の特性を十分考慮する必要がある。
【０００４】
例えば、非特許文献１には、高分子電解質膜燃料電池のカソードについて、２次元的、非等温、気相液相二相、複数要素モデルを用いた電池の電圧電流解析が述べられている。ここでは、気体と液体の流れについてダルシー（Ｄａｒｃｙ）の式を用い、接触角、透過率、多孔率等を実験で求め、毛管圧と飽和度の関係は水銀ポロジメトリに基づいた方法を適用し、相対透過率は、このデータからバーダイン（Ｂｕｒｄｉｎｅ）の式を適用して求めることが述べられている。
【０００５】
【非特許文献１】Ｍ．アコスタ等，高分子電解質膜燃料電池のカソードにおける２次元的、非等温、二相、複数要素流れのモデル化，ジャーナル  オブ  パワーソース，２００６年，１５９号，ｐ１１２３−１１４１
【発明の開示】
【発明が解決しようとする課題】
【０００６】
非特許文献１の方法によれば、ダルシーの式において必要な相対透過率を、その多孔体についての実験結果とバーダインの式を用いて算出している。これによって、その多孔体の圧力損失をシミュレーションで求めることができる。
【０００７】
非特許文献１における実験結果では、多孔体の気孔率分布、すなわち気孔径の頻度分布は、必ずしも１つの気孔径に集まっているわけではなく、ばらつきがあるが、これを適当な分布形を有するものとして取り扱っている。
【０００８】
実際の多孔体を調べてみると、非特許文献１で述べられているような気孔率分布とは別に、気孔率分布が複数のピークを有するものがある。複数のピークを有する場合でも、その分布形を用いて非特許文献１の方法を適用することができる。これを実際のシミュレーションで実行してみると、圧力損失の計算結果が実験結果と合わないことがある。
【０００９】
本発明の目的は、多孔体の気孔率分布が複数の気孔率ピークを有するときに、実験結果との適合性のよい計算結果を出すことを可能とする気液二相流れ評価シミュレーション装置及び気液二相流れ評価シミュレーションプログラムを提供することである。
【課題を解決するための手段】
【００１０】
本発明に係る気液二相流れ評価シミュレーション装置は、多孔体における気液二相流れについて、多孔体の気孔率と、多孔体の構造で定まる絶対透過率Ｋと、多孔体における気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ与えて、所定の関係式に基づいて、多孔体の流入側と流出側との間の差圧である圧力損失ΔＰを求める気液二相流れ評価シミュレーション装置であって、多孔体の気孔率分布が複数の気孔率ピークを有するとき、多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとしてモデル化する手段と、各層ごとに、その気孔率分布に応じて、絶対透過率Ｋと、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ割り当てる手段と、積層モデルの各層のうち、気液二相流れが多孔体に供給される流入側層を第１層とし、流入圧力ＰINを与えて、第１層の圧力損失ΔＰ1から第１層の流出圧力側のＰOUT1を求める手段と、積層モデルの各層のうち、第１層の次に流出側に配置された層を第２層とし、第１層の流出圧力ＰOUT1を第２層に対する流入圧力ＰIN2として与え、第２層の圧力損失ΔＰ2から第２層の流出圧力ＰOUT2を求める手段と、積層モデルを構成する層数に応じてこれを繰り返し、積層モデルの各層のうち、気液二相流れが多孔体から流出する流出側層について求められた流出圧力をＰOUTとし、ΔＰ＝ＰIN−ＰOUTを多孔体の圧力損失ΔＰとして出力する手段と、を備えることを特徴とする。
【００１１】
また、本発明に係る気液二相流れ評価シミュレーション装置は、気体密度ρGと気体粘性係数μGと気体質量速度ＧGを有する気体と、液体密度ρLと気体粘性係数μGと液体質量速度ＧLを有する液体とが気液二相流れとして、厚みがＬである多孔体に供給されるとき、入力条件として、多孔体の気孔率と、多孔体の構造で定まる絶対透過率Ｋと、を与え、気孔率に基づいて液体の飽和度ｓと毛管圧ｐCを求め、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式に基づいて、気体と液体のそれぞれについて、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ求め、ΔＰ／Ｌ＝（μGＧG）／（ｋGＫρG）＝（μLＧL）／（ｋLＫρL）の関係式に基づき圧力損失ΔＰを求める気液二相流れ評価シミュレーション装置であって、多孔体の気孔率分布が複数の気孔率ピークを有するとき、多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとし、各層をｉ層として、ｉ層のそれぞれについて、その気孔率分布に応じて、厚さＬiと絶対透過率Ｋiとを割り当てるモデル化手段と、ｉ層のそれぞれについて、その気孔率分布に応じて、気体相対透過率ｋGiと液体相対透過率ｋLiとを求める相対透過率算出手段と、積層モデルの各層のうち、気液二相流れが多孔体に供給される流入側層を第１層とし、流入圧力ＰINを与え、ｉ＝１として、ΔＰ／Ｌ＝（ＰIN−ＰOUT1）／Ｌ1＝（μGＧG）／（ｋG1Ｋ1ρG）＝（μLＧL）／（ｋL1Ｋ1ρL）から、第１層の流出圧力ＰOUT1を求める手段と、積層モデルの各層のうち、第１層の次に流出側に配置された層を第２層とし、第１層の流出圧力ＰOUT1を第２層に対する流入圧力ＰIN2として与える手段と、第２層について、流入圧力ＰIN2＝ＰOUT1を与え、ｉ＝２として、ΔＰ／Ｌ＝（ＰIN2−ＰOUT2）／Ｌ2＝（μGＧG）／（ｋG2Ｋ2ρG）＝（μLＧL）／（ｋL2Ｋ2ρL）から、第２層の流出圧力ＰOUT2を求める手段と、積層モデルを構成する層数に応じてこれを繰り返し、積層モデルの各層のうち、気液二相流れが多孔体から流出する流出側層について求められた流出圧力をＰOUTとし、ΔＰ＝ＰIN−ＰOUTを多孔体の圧力損失ΔＰとして出力する手段と、を備えることを特徴とする。
【００１２】
また、本発明に係る気液二相流れ評価シミュレーション装置において、相対透過率算出手段は、各ｉ層について、ｉ層の気孔率分布を構成する各気孔径ＤPについて液体の飽和度ｓiを割り当てて飽和度分布ｓi＝ｆi（Ｄp）を求め、この飽和度分布ｓiに対応する毛管圧分布ｐCi＝ｇi（ｓ）を求める手段と、各ｉ層についての飽和度分布ｓi＝ｆi（Ｄp）と毛管圧分布ｐCi＝ｇi（ｓ）とに基づき、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式から、各ｉ層についての気体相対透過率の分布ｋG＝ｋG（ｓ）と液体相対透過率の分布ｋL＝ｋL（ｓ）を求める手段と、を含むことが好ましい。
【００１３】
また、本発明に係る気液二相流れ評価シミュレーション装置において、ｉ層を複数の分割要素に分割し、ｉ層の気体相対透過率の分布ｋG＝ｋG（ｓ）と液体相対透過率の分布ｋL＝ｋL（ｓ）に基づいて、各分割要素を分割要素Ｊとして、分割要素Ｊのそれぞれに気体相対透過率ｋGJと液体相対透過率ｋLJを割り当て、有限要素法を用いて多孔体の圧力損失ΔＰを求めることが好ましい。
【００１４】
また、本発明に係る気液二相流れ評価シミュレーションプログラムは、気体密度ρGと気体粘性係数μGと気体質量速度ＧGを有する気体と、液体密度ρLと気体粘性係数μGと液体質量速度ＧLを有する液体とが気液二相流れとして、厚みがＬである多孔体に供給されるとき、入力条件として、多孔体の気孔率と、多孔体の構造で定まる絶対透過率Ｋと、を与え、気孔率に基づいて液体の飽和度ｓと毛管圧ｐCを求め、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式に基づいて、気体と液体のそれぞれについて、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ求め、ΔＰ／Ｌ＝（μGＧG）／（ｋGＫρG）＝（μLＧL）／（ｋLＫρL）の関係式に基づき圧力損失ΔＰを求める気液二相流れ評価シミュレーション装置において実行される気液二相流れ評価シミュレーションプログラムであって、多孔体の気孔率分布が複数の気孔率ピークを有するとき、多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとし、各層をｉ層として、ｉ層のそれぞれについて、その気孔率分布に応じて、厚さＬiと絶対透過率Ｋiとを割り当てるモデル化処理手順と、ｉ層のそれぞれについて、その気孔率分布に応じて、気体相対透過率ｋGiと液体相対透過率ｋLiとを求める相対透過率算出処理手順と、積層モデルの各層のうち、気液二相流れが多孔体に供給される流入側層を第１層とし、流入圧力ＰINを与え、ｉ＝１として、ΔＰ／Ｌ＝（ＰIN−ＰOUT1）／Ｌ1＝（μGＧG）／（ｋG1Ｋ1ρG）＝（μLＧL）／（ｋL1Ｋ1ρL）から、第１層の流出圧力ＰOUT1を求める処理手順と、積層モデルの各層のうち、第１層の次に流出側に配置された層を第２層とし、第１層の流出圧力ＰOUT1を第２層に対する流入圧力ＰIN2として与える処理手順と、第２層について、流入圧力ＰIN2＝ＰOUT1を与え、ｉ＝２として、ΔＰ／Ｌ＝（ＰIN2−ＰOUT2）／Ｌ2＝（μGＧG）／（ｋG2Ｋ2ρG）＝（μLＧL）／（ｋL2Ｋ2ρL）から、第２層の流出圧力ＰOUT2を求める処理手順と、積層モデルを構成する層数に応じてこれを繰り返し、積層モデルの各層のうち、気液二相流れが多孔体から流出する流出側層について求められた流出圧力をＰOUTとし、ΔＰ＝ＰIN−ＰOUTを多孔体の圧力損失ΔＰとして出力する処理手順と、が実行される。
【発明の効果】
【００１５】
上記構成の少なくとも１つにより、気液二相流れ評価シミュレーション装置は、多孔体における気液二相流れについて、多孔体の気孔率と、多孔体の構造で定まる絶対透過率Ｋと、多孔体における気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ与えて、所定の関係式に基づいて、多孔体の流入側と流出側との間の差圧である圧力損失ΔＰを求める気液二相流れ評価シミュレーション装置である。この基本構成は、多孔体の気孔率分布が複数の気孔率ピークを有するときの配慮を除けば、非特許文献１に開示されるものと同じである。
【００１６】
上記構成においては、多孔体の気孔率分布が複数の気孔率ピークを有するときに、多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとしてモデル化し、各層ごとに、その気孔率分布に応じて、絶対透過率Ｋと、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ割り当てる。
【００１７】
そして、気液二相流れが多孔体に供給される流入側層を第１層とし、流入圧力ＰINを与えて、第１層に割り当てられたデータに基づき、上記基本構成の手順に従って、この第１層の流出圧力を求める。
【００１８】
第２層については、第１層の流出圧力を第２層の流入圧力として、同様に、第２層に割り当てられたデータに基づき、上記基本構成の手順に従って、この第２層の流出圧力を求める。すなわち、第１層の流出圧力＝第２層流入圧力が境界条件となる。
【００１９】
これをモデル化された複数の層について繰り返し、最後の層である流出層について、その前の層の流出圧力を流出層の流入圧力として、同様に、流出層に割り当てられたデータに基づき、上記基本構成の手順に従って、この流出層の流出圧力ＰOUTを求める。
【００２０】
このようにして、モデル化された多孔体の流出圧力ＰOUTが求められるので、圧力損失ΔＰ＝ＰIN−ＰOUTを求めることができる。求められた計算結果は、積層モデル化しない場合に比べ、実験結果とよく一致することが確かめられた。したがって、上記手順を実行することで、多孔体の気孔率分布が複数の気孔率ピークを有するときに、実験結果と適合性のよい気液二相流れ評価シミュレーションを行うことができる。
【００２１】
また、絶対透過率Ｋと相対透過率ｋ等を与えて多孔体の圧力損失を求める式として、ΔＰ／Ｌ＝（μGＧG）／（ｋGＫρG）＝（μLＧL）／（ｋLＫρL）の関係式を用いることができる。この式は、流体の運動量保存則から導かれる関係式で、ダルシーの関係式と呼ばれる。また、相対透過率を液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式に基づいて求めることができる。この関係式としては、バーダインの式が知られている。このように、既に実績のあるダルシーの関係式、バーダインの式を用いることで、多孔体の圧力損失について実験結果と適合性のよい気液二相流れ評価シミュレーションを行うことができる。
【００２２】
また、モデル化した各層について、気孔率分布を構成する各気孔径ＤPについて液体の飽和度ｓiを割り当てて飽和度分布ｓi＝ｆi（Ｄp）を求める。また、この飽和度分布ｓiに対応する毛管圧分布ｐCi＝ｇi（ｓ）を求める。関数ｆi，ｇiの関数形は、予め実験等で求めることができるので、各層の気孔径ＤPの分布を与えることで、各層のｓi，ｐCiを求めることができる。そして、バーダインの式として知られている液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係式から、各ｉ層についての気体相対透過率の分布ｋG＝ｋG（ｓ）と液体相対透過率の分布ｋL＝ｋL（ｓ）を求める。このように、各層の気孔径分布に応じて、各層の相対飽和度を求めるので、実験結果と適合性のよい気液二相流れ評価シミュレーションを実行することができる。
【００２３】
また、モデル化された各ｉ層を複数の分割要素に分割し、ｉ層の気体相対透過率の分布ｋG＝ｋG（ｓ）と液体相対透過率の分布ｋL＝ｋL（ｓ）に基づいて、有限要素法を用いて多孔体の圧力損失ΔＰを求める。したがって、きめ細かな気液二相流れ評価シミュレーションを行うことができる。
【発明を実施するための最良の形態】
【００２４】
以下に図面を用いて、本発明に係る実施の形態につき、詳細に説明する。以下では、多孔体の気孔率分布が２つのピークを有するものとして説明するが、これは説明のための一例であって、３以上のピークを有するものであってもよい。シミュレーションの対象としての多孔体は、気孔を有するポーラスな材料であればよく、多孔質フィルム、多孔質膜あるいは多孔質の焼結材等のほか、半透過膜、中空糸等であってもよい。また、以下では、気液二相流れとして、空気と水から構成される流体を説明するが、気体は空気以外のガスであってもよく、液体は水以外のものであってもよい。また、以下における数値等は説明のための一例であって、多孔体の種類、供給される流体の種類等に応じ、適宜変更が可能である。
【００２５】
図１は、気液二相流れ評価シミュレーション装置１０の構成を示す図である。気液二相流れ評価シミュレーション装置１０は、多孔体に気液二層流れを与えたときに、多孔体と気液二相流れの条件を与えて、多孔体における圧力損失ΔＰをシミュレーションで求める装置である。ここでは特に、多孔体の気孔率分布が複数のピークを有するときの圧力損失を求めるシミュレーションを実行する機能を有する。
【００２６】
ここでまず、多孔体に気液二相流れが供給されたときの多孔体の両端における圧力損失を求める方法について、一般的な説明を行い、その後に、気液二相流れ評価シミュレーション装置１０の構成を説明する。
【００２７】
多孔体とは、一般的に、気孔を有し、気体と液体を透過する性質を有する物質である。また、気液二相流れとは、気体と液体とから構成される流れのことで、例えば、気体として空気を、液体として水を構成要素とする流れである。このような流体を、気液二相流体と呼んでもよい。
【００２８】
多孔体と気液二相流れの条件を与えて、多孔体における圧力損失ΔＰを求めるときに基礎となる関係式は、ダルシーの関係式、あるいはダルシー則と呼ばれるものである。この関係式は、流体の運動量保存則等から導かれるもので、流速ｕと粘性係数μを有する流体が、絶対透過率Ｋ、相対透過率ｋを有する物質をｄｘの距離を通過すると、圧力損失ｄＰは、ｄＰ＝［（μ・ｕ）／Ｋ・ｋ］ｄｘで与えられるというものである。ここで、絶対透過率Ｋは、多孔体等の物質の構造で定まる特性係数で、次元は［（面積）2］で、例えば［ｍ2］である。相対透過率ｋは、絶対透過率Ｋに対する割合を示す係数で、０から１の値をとる。例えば、流体が水の場合には、多孔体の含水率に対応するものである。
【００２９】
図２に、気液二相流れの場合のダルシーの関係式を示す。ここでは、厚みがｄｘである多孔体３０に、気液二相流れ３２として、気体流速ｕGと気体粘性係数μGを有する気体（Ｇ）と、液体流速ｕLと液体粘性係数μLを有する液体（Ｌ）とで構成されるときに、多孔体３０の両端における圧力差ｄＰを表す式が、気体（Ｇ）と液体（Ｌ）のそれぞれについて示されている。
【００３０】
なお、図２に示されるように、流速ｕは、ｕ＝Ｇ／ρで表すこともできる。ここで、ρは密度で、Ｇは質量速度である。質量速度とは、［質量／（面積×時間）］の次元を有し、流れにおける単位面積当り、単位時間当たりの質量である。図２では、気体質量速度ＧG、液体質量速度ＧL、気体密度ρG、液体密度ρLと、気体流速ｕG、液体流速ｕLとの関係が示されている。
【００３１】
したがって、多孔体３０について、厚みｄｘ、絶対透過率Ｋ、気体の相対透過率ｋG、液体の相対透過率ｋLを与え、気液二相流れ３２について気体流速ｕG、気体粘性係数μG、液体流速ｕL、液体粘性係数μLを与えることで、多孔体３０の両端における圧力損失ｄＰを求めることができる。なお、ここでは、気体（Ｇ）における圧力損失（差圧）と液体における圧力損失（差圧）を等しいものとしている。
【００３２】
上記の各パラメータが、一定値であるときは、上記の式により、多孔体の両端における
圧力損失を容易に求めることができる。各パラメータが、分布を有するときでも、分布形が比較的簡単な関数形で与えられるときには、例えば、有限要素法等を用いてシミュレーションを進めることができる。また、分布形が複雑な関数形であっても、近似が可能なときは、近似的な関数形を用いてシミュレーションを進めることができる。
【００３３】
図２の関係式において、多孔体の構造等に関係して複雑な関数形となるのが相対透過率ｋである。気液二相流れにおいては、気体の相対透過率ｋG、液体の相対透過率ｋLの２つが必要であり、これらは相互に関連がある。相対透過率ｋについては、バーダインの関係式が知られている。バーダインの関係式は、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの間を関係づけるもので、たとえば、上記の非特許文献１に示されている。なお、毛管圧ｐCは、飽和度ｓの関数であることが知られている。
【００３４】
したがって、多孔体について、液体の毛管圧ｐCと飽和度ｓの関数形と、飽和度ｓが分かれば、気体の相対透過率ｋG、液体の相対透過率ｋLが分かる。ここで、飽和度ｓが分布形を有するときは、その分布形について、バーダインの関係式を適用することになる。飽和度ｓは、上記のように含水率に対応するものである。したがって、多孔体が一様で均一な気孔からなる場合には飽和度ｓは一定値となるが、多孔体が不均一な気孔径分布を有するときは、飽和度ｓも分布形を有することになり、バーダインの関係式の適用も複雑になる。
【００３５】
したがって、多孔体の気孔率分布の関数形が複雑であればあるほど、バーダインの関係式の計算が複雑になり、気体の相対透過率ｋG、液体の相対透過率ｋLを求めることが困難になる。そこで、気孔率分布の関数形を近似式で表現することが考えられるが、近似式の内容によっては、気体の相対透過率ｋG、液体の相対透過率ｋLが実際のものと乖離することが生じる。例えば、多孔体の気孔率分布が複数のピークを有する場合、平均的な気孔率分布として近似すると、計算で求められる気体の相対透過率ｋG、液体の相対透過率ｋLが実際のものと乖離し、これに基づく多孔体の圧力損失の計算結果も実験結果と乖離する。これが本発明に関連する課題である。
【００３６】
ここで再び図１に戻り、気液二相流れ評価シミュレーション装置１０の構成を説明する。気液二相流れ評価シミュレーション装置１０（以下では、簡単にシミュレーション装置１０と呼ぶことにする。）は、演算処理等を実行するＣＰＵ１２と、計算のための条件を入力するキーボード等の入力部１４と、計算結果を出力するディスプレイ等の出力部１６と、プログラム等を記憶する記憶装置１８とを有して構成される。なお、これらの要素が互いに内部バスで接続される。かかる気液二相流れ評価シミュレーション装置１０は、有限要素法の実行を含むシミュレーション計算に適したコンピュータで構成することができる。
【００３７】
ＣＰＵ１２は、多孔体の気孔率分布が複数の気孔率ピークを有するとき、多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとしてモデル化するモデル化モジュール２０を有する。ここでは、多孔体の気孔率分布が２つの気孔率ピークを有するとして、これを２つの層が積層されたモデルとし、気液二相流れが供給される側を流入層とし、気液二相流れが排出される側を流出層とする。
【００３８】
また、ＣＰＵ１２は、モデル化された各層について、気体相対透過率ｋGと液体相対透過率ｋLとをそれぞれ算出する相対透過率算出モジュール２２を有する。なお、相対透過率算出は、各層の気孔率分布に基づいて、バーダインの関係式によって実行される。求められた気体相対透過率ｋGと液体相対透過率ｋLは、対応する各層に割り当てられる。なお、各層には、その気孔率分布に応じた絶対透過率Ｋも割り当てられる。
【００３９】
また、ＣＰＵ１２は、流入層について、流入圧力ＰINを与えたときの流入層における流出圧力ＰOUT1を算出する流入層ＰOUT1算出モジュール２４を有する。ここでは、流入層に割り当てられた絶対透過率Ｋと、気体相対透過率ｋGと液体相対透過率ｋLとに基づき、ダルシーの関係式を用いて、流入層の両端の圧力損失ΔＰを求め、ΔＰ＝ＰIN−ＰOUT1から、流入層の流出圧力のＰOUT1が算出される。
【００４０】
また、ＣＰＵ１２は、多孔体の流出圧力ＰOUTとなる流出層の流出圧力を算出する流出層ＰOUT算出モジュール２６を有する。ここでは、流入層と流出層の境界条件として、流入層ＰOUT1を流出層の流入圧力とする。そして、流出層に割り当てられた絶対透過率Ｋと、気体相対透過率ｋGと液体相対透過率ｋLとに基づき、ダルシーの関係式を用いて、流出層の両端の圧力損失ΔＰを求め、ΔＰ＝ＰOUT1−ＰOUTから、流出層の流出圧力、すなわち多孔体の流出圧力ＰOUTが算出される。
【００４１】
そして、ＣＰＵ１２は、多孔体の両端の圧力損失ΔＰを求める圧力損失ΔＰ算出モジュール２８を有する。ここでは、ΔＰ＝ＰIN−ＰOUTが演算されて出力される。
【００４２】
ここで、上記各機能の実行に必要な入力条件等は、入力部１４からの入力を受け取って取得され、ダルシーの関係式、バーダインの関係式等を含む演算処理手順を与えるプログラム等は、記憶装置１８に格納されて、必要に応じ読み出される。
【００４３】
かかる各機能は、ソフトウェアで実現でき、具体的には、対応する気液二相流れ評価シミュレーションプログラムを実行することで実現される。
【００４４】
次に、上記構成の作用を図３、図７のフローチャート、及び図４から図６、図８から図１０の各説明図を用いて説明する。図３は、気孔率分布が複数のピークを有する多孔体に、気液二相流れが供給されたときの多孔体の両端の圧力損失を算出する手順を示すフローチャートで、各手順は、気液二相流れ評価シミュレーションプログラム評価シミュレーションプログラムの各処理手順にそれぞれ対応する。また、図７は、図３における相対透過率算出工程（Ｓ１２）の内部手順を示すフローチャートである。図４から図６はモデル化を説明する図、図８から図１０は、相対透過率の算出の様子を説明する図である。
【００４５】
気孔率分布が複数のピークを有する多孔体に、気液二相流れが供給されたときの多孔体の両端の圧力損失を算出するには、まず多孔体を複数の層が積層されたものとするモデル化が行われる。図３では、一般的に複数層として説明するために、ｉ層モデル化として示されている（Ｓ１０）。この工程は、ＣＰＵ１２のモデル化モジュール２０の機能によって実行される。モデル化は、多孔体の気孔率分布が複数の気孔率ピークを有するとき、多孔体を、各気孔率ピークに対応する気孔率分布を有する各層が順次積層配置されて構成される積層モデルとし、各層をｉ層として、ｉ層のそれぞれについて、その気孔率分布に応じて、厚さＬiと絶対透過率Ｋiとを割り当てることで行われる。
【００４６】
モデル化の様子を図４から図６を用いて説明する。図４は、多孔体の気孔径分布の様子を示す一例で、左側の図は、横軸に気孔径ＤPをとり、縦軸に気孔率の累計をとったものである。なお、気孔率の累計は、気孔径ＤPの大きい方からはじめ、気孔率ＤPの小さい方に向かって累計が行われる。この気孔率累計の方法は、水銀ポロジメトリ法を用いるもので、水銀ポロジメトリ法によれば、気孔径ＤPの大きい方から水銀が満たされてゆくからである。
【００４７】
図４の左側の図では、気孔率累計曲線が２箇所の気孔径のところで変化が大きくなっている。ここでは、Ａ，Ｂでこの２箇所が示されている。これは、この２箇所Ａ，Ｂの気孔径に対応する気孔率が高いことを示している。このことは、気孔率累計曲線を気孔径ＤPで微分した様子を示す図４の右側の図で見るとよく分かる。気孔率累計曲線を気孔径ＤPで微分したものは、気孔率分布となるので、図４の右側の図で示されるように、この多孔体の場合は、気孔率分布が２つのピークを有する。ここでも、２つのピークに対応する箇所をＡ，Ｂで示されている。すなわち、この多孔体の場合、Ａの箇所、すなわち、小さい気孔径のところでピークを有する第１の気孔率分布と、Ｂの箇所、すなわち、大きい気孔径のところでピークを有する第２の気孔率分布を有していることになる。
【００４８】
図５は、気孔率分布が２つのピークを有する多孔体を模式的に示す図である。このように、この多孔体は、大きい気孔径を有する大気孔径グループと、小さい気孔径を有する小気孔径グループとを有する。気孔率分布が２つのピークを有する典型的な場合は、大気孔径グループに属する気孔のうちで最も小さい気孔径は、小気孔径グループに属する気孔のうちで最も大きな気孔径よりも、さらに大きい気孔径である場合である。勿論、いくらかの重なりがあってもよいが、少なくとも、図４の右側の図に示されるように、気孔率分布において、ピークが２つある場合である。
【００４９】
このように、気孔率分布が２つのピークを有する場合には、この多孔体を小気孔径グループに属する気孔で構成される第１層と、大気孔径グループに属する気孔で構成される第２層の２つの層が積層されたものとしてモデル化する。その様子を図６に示す。ここでは、ｉ＝１として、小気孔径グループに属する気孔で構成される第１層が、ｉ＝２として、大気孔径グループに属する気孔で構成される第２層がそれぞれ示される。
【００５０】
なお、気孔率分布が３以上の場合には、そのピーク数に応じた層数の積層としてモデル化が行われる。以下では、特に断らない限り、気孔率分布が２つのピークを有する場合について説明を続ける。
【００５１】
モデル化に際して、もともとの多孔体が厚みＬと絶対透過率Ｋを有しているとして、モデル化された第１層と第２層に、それぞれの厚みＬ1，Ｌ2と、それぞれの絶対透過率Ｋ1，Ｋ2が割り振られる。厚みについては、Ｌ＝Ｌ1＋Ｌ2として、Ｌ1，Ｌ2との間では、小気孔径グループに属する気孔の数と、大気孔径グループに属する気孔の数とを比較して、割り振りが行われる。すなわち、小気孔径グループの気孔率と大気孔径グループの気孔率とに応じて、それぞれの厚みＬ1，Ｌ2が割り振られる。絶対透過率については、小気孔径の気孔で構成されるときの絶対透過率と、大気孔径の気孔で構成されるときの絶対透過率が異なるので、それぞれの値が割り振られる。すなわち、小気孔径グループのピークの気孔径と、大気孔径グループのピークの気孔径に応じて、それぞれの絶対透過率Ｋ1，Ｋ2が割り振られる。
【００５２】
再び図３に戻り、次に相対透過率が算出される（Ｓ１２）。この工程は、ＣＰＵ１２の相対透過率算出モジュール２２の機能によって実行される。相対透過率算出工程は、図７に示されるサブ工程に分けることができる。
【００５３】
図７において、最初は、気孔率分布から飽和度ｓ＝ｆ（ＤP）を求める工程である。気孔率分布とは、図４の右側の分布のことで、例えば、小気孔径グループの場合、ピークとなる気孔径を中心として、気孔径にばらつきがある。大気孔径グループの場合も同様である。このように、同じ気孔径グループの中での気孔径の分布のことを、ここで気孔率分布と称している。
【００５４】
図８は、小気孔径グループにおける気孔率累計曲線と、大気孔径グループにおける気孔率累計曲線を示す図である。図８（ａ）には、小気孔径グループにおける気孔率累計曲線が実線で、もともとの多孔体全体の気孔率累計曲線が破線で示されている。同様に、図８（ｂ）には、大気孔径グループにおける気孔率累計曲線が実線で、もともとの多孔体全体の気孔率累計曲線が破線で示されている。なお、もともとの多孔体全体の気孔率累計曲線
は、図４の左側の気孔率累計曲線と同じである。そして、小気孔径グループの気孔率累計曲線は、もともとの多孔体全体の気孔率累計曲線の中で、小気孔径グループの気孔径に相当する部分を抜き出したものであり、大気孔径グループの気孔率累計曲線は、もともとの多孔体全体の気孔率累計曲線の中で、大気孔径グループの気孔径に相当する部分を抜き出したものである。
【００５５】
ここで、上記のように、気孔率累計曲線は、水銀ポロジメトリ法によって求められている。すなわち、気孔率累計曲線は、気孔が水銀で飽和してゆく特性を示している。このことから、図８の気孔率累計曲線は、水銀の飽和度と気孔径ＤPとの関係を示していることになる。そこで、この気孔率累計曲線の関数形を、液体の飽和度ｓと気孔径ＤPとの関係を示す関数形であるとする。
【００５６】
具体的には、小気孔径グループの場合、図８（ａ）の実線の曲線の最大値をｓ＝１に規格化した関数形を求め、これを小気孔径グループの層、すなわちｉ＝１の層における飽和度ｓの気孔径ＤPについての関数形とする。図８（ａ）ではｓ1＝ｆ1（ＤP）として示されている。上記のように、この関数形は、ＤPが図８（ａ）の横軸の一番左側のところでｓ＝１となるように規格化されている。
【００５７】
同様に、大気孔径グループの場合、図８（ｂ）の実線の曲線の最大値をｓ＝１に規格化した関数形を求め、これを大気孔径グループの層、すなわちｉ＝２の層における飽和度ｓの気孔径ＤPについての関数形とする。図８（ｂ）ではｓ2＝ｆ2（ＤP）として示されている。上記のように、この関数形は、ＤPが図８（ａ）の横軸の一番左側のところでｓ＝１となるように規格化されている。
【００５８】
再び図７に戻り、各層について気孔径ＤPに対する飽和度ｓの関数形が求められると、この関数形を用いて、各層について、毛管圧ｐCと飽和度ｓの関数形を求める（Ｓ２４）。この関数形を求めるには、表面張力σと接触角θと毛管圧ｐCと気孔径ＤPの関係式であるＤP＝４σｃｏｓθ／ｐCを用いる。その様子が、図９に示される。図９（ａ）には、ｉ＝１、すなわち小気孔径グループにおける毛管圧ｐC1と飽和度ｓの関数形として、ｐC1＝ｇ1（ｓ）が示されている。同様に、図９（ｂ）には、ｉ＝２、すなわち大気孔径グループにおける毛管圧ｐC2と飽和度ｓの関数形として、ｐC1＝ｇ1（ｓ）が示されている。
【００５９】
再び図７に戻り、次に、バーダインの関係式にもとづき、各層について、気体相対透過率ｋGと液体相対透過率ｋLを飽和度ｓの関数形として求める（Ｓ２６）。バーダインの関係式は、上記のように、液体の飽和度ｓと毛管圧ｐCと相対透過率ｋとの関係を示すものである。したがって、各層について、Ｓ２４において求められた毛管圧ｐCiと飽和度ｓとの間の関数形ｇiを用いて、各層について気体相対透過率ｋGと液体相対透過率ｋLを飽和度ｓの関数として求めることができる。
【００６０】
その様子を図１０に示す。図１０（ａ）は、ｉ＝１、すなわち小気孔径グループにおける気体相対透過率ｋGと液体相対透過率ｋLが飽和度ｓの関数としてそれぞれ示されている。同様に、図１０（ｂ）には、ｉ＝２、すなわち大気孔径グループにおける気体相対透過率ｋGと液体相対透過率ｋLが飽和度ｓの関数としてそれぞれ示されている。このようにして、各層について、気体相対透過率ｋGと液体相対透過率ｋLが飽和度ｓの関数として求められる。
【００６１】
ここで再び図３に戻る。上記のように、各層について、気体相対透過率ｋGと液体相対透過率ｋLが飽和度ｓの関数として求められると、各層に、気体相対透過率ｋGと液体相対透過率ｋLが飽和度ｓの関数として割り当てられる。上記のように飽和度ｓは、気孔径ＤPの関数であるので、気孔径ＤPが定まれば、飽和度ｓが定まり、飽和度ｓが定まれば、気体相対透過率ｋGと液体相対透過率ｋLが定まる。
【００６２】
このようにして、各層について気体相対透過率ｋGと液体相対透過率ｋLが定まると、各層の圧力損失がダルシーの関係式を用いて算出される。最初は、ｉ＝１の層について、流入圧力ＰINが与えられて、その圧力損失が求められ、その層の流出圧力ＰOUT1が算出される（Ｓ１４）。この工程は、ＣＰＵ１２の流入層ＰOUT1算出モジュール２４の機能によって実行される。上記の例では、２層のモデル化を行っているので、ｉ＝１の層、すなわち、小気孔径グループの第１層について、その層の流出圧力ＰOUT1が算出される。
【００６３】
ダルシーの関係式は、質量速度Ｇと密度ρを用いて表すと、次のような式である。すなわち、多孔体の厚みをＬ、圧力損失をΔＰとし、絶対透過率Ｋ、気体密度ρG、気体粘性係数μG、気体相対透過率ｋG、液体密度ρL、液体粘性係数μL、液体相対透過率ｋLを与えると、ΔＰ／Ｌ＝（μGＧG）／（ｋGＫρG）＝（μLＧL）／（ｋLＫρL）となる。ここで、気体密度ρG、気体粘性係数μG、液体密度ρL、液体粘性係数μLは、気液二相流れの特性値として、図１の入力部１４から与えられる。
【００６４】
したがって、ｉ＝１の層について、厚みＬ1、絶対透過率Ｋ1はＳ１０で述べたようにして割り当てられ、気体相対透過率ｋG、液体相対透過率ｋLは上記のようにＳ１２において算出されて割り当てられるので、ダルシーの関係式に必要なパラメータが全部揃ったことになる。したがって、ダルシーの関係式に基づいて圧力損失ΔＰを求めることができる。具体的には、流入圧力ＰINが与えられて、ダルシーの関係式からΔＰが求められると、ΔＰ＝ＰIN−ＰOUT1であるので、その層の流出圧力ＰOUT1が算出される。
【００６５】
次にｉ＝２の層について、同様の計算が実行される。上記の例では、２層のモデル化を行っているので、ｉ＝２の層、すなわち、大気孔径グループの第２層について同様の計算が実行される。この場合、ｉ＝２層に対する流入圧力は、ｉ＝１の層の流出圧力ＰOUT1が用いられる。すなわち、ｉ＝１の層とｉ＝２の層とは積層されているモデルであるので、境界条件として、ある層の流出圧力は、隣接する層の流入圧力であるとするものである。したがって、ｉ＝２の層については、流入圧力＝ＰOUT1として、その層の流出圧力が求められる（Ｓ１６）。
【００６６】
そして、上記の例では、２層のモデル化を行っているが、３層以上のモデル化の場合には、ｉ＝２の層の流出圧力をｉ＝３の層の流入圧力とされる。以後、このような境界条件の下で、モデル化された各層について順次その流出圧力が求められる。最後の層について求められた流出圧力は、モデル化の元の多孔体にとっての流出圧力であるので、これがＰOUTとして算出される（Ｓ１８）。Ｓ１６からＳ１８の工程は、ＣＰＵ１２の流出層ＰOUT算出モジュール２６の機能によって実行される。上記の例では、２層のモデル化であるので、Ｓ１６の内容がＳ１８となり、ｉ＝２の層の流出圧力がＰOUTとなる。
【００６７】
そして、Ｓ１４におけるＰINとＳ１８において算出されたＰOUTから、モデル化の元の多孔体の両端における圧力損失ΔＰがΔＰ＝ＰIN−ＰOUTとして求められ、出力される（Ｓ２０）。この工程は、ＣＰＵ１２の圧力損失ΔＰ算出モジュール２８の機能によって実行される。
【００６８】
上記では、気孔径ＤPが定まれば、飽和度ｓが定まり、飽和度ｓが定まれば、気体相対透過率ｋGと液体相対透過率ｋLが定まるので、各層の圧力損失を求められるとして説明した。実際には、図４に示されるように、小気孔径グループにおいて気孔径ＤPは１つの値ではなく、ピークの気孔径を中心として分布を有する。大気孔径グループにおいても同様である。多孔体を複数の層が積層されるものとしてモデル化したときに、各層について気孔径がある分布を有する場合には、例えば有限要素法を用いて、上記と同様なシミュレーションを実行することができる。
【００６９】
すなわち、モデル化した各層を複数の分割要素に分割し、各分割要素には１つの気孔径、あるいは飽和度を割り当てるようにする。気孔径ＤPあるいは飽和度ｓの割り当ては、図４で求められた気孔率分布あるいは図８で説明した飽和度分布に基づいて行うことができる。このようにして、各分割要素に１つの気孔径ＤPあるいは１つの飽和度ｓが割り当てられると、各分割要素にはダルシーの関係式に必要なパラメータが全部揃ったことになるので、Ｓ１４等で説明したように、流入圧力を与えれば、その分割要素について流出圧力を求めることができる。
【００７０】
そして、各層について、その層の流入側から流出側に向かって、Ｓ１４からＳ１８で説明したと同様の処理を順次実行する。すなわち、流入側にある最初の分割要素にはその層の初期条件である流出圧力を与え、その分割要素について求められた流出圧力を次の分割要素の流入圧力とし、その分割要素で求められた流出圧力をさらに次の分割要素の流入圧力とする。これを順次繰り返し、最後の流出側の分割要素で求められた流出圧力が、その層の流出圧力となる。
【００７１】
このように、各層が気孔率のピークを中心に気孔率分布を有する場合でも、分割要素の大きさを気孔率分布に合わせて選択することで、有限要素法等を用いて、各層の圧力損失を求めることができる。
【００７２】
図１１は、複数の層にモデル化してＳ１２で説明した方法で求めた気体相対透過率ｋGと液体相対透過率ｋLの計算結果と、別途実験で求めた値とを示す図である。計算結果は実線で示され、実験結果は○印のデータ点で示される。なお、破線は、２つのピークを有する気孔率分布を１つの平均的気孔率分布で近似した場合の計算結果である。図１１に示されるように、気孔率分布が複数のピークを有する場合に、複数の層にモデル化して計算を行うことで、気体相対透過率ｋGと液体相対透過率ｋLの計算結果が実験値とよく一致することが分かる。
【００７３】
図１２と図１３は、多孔体の圧力損失ΔＰと、液体（Ｌ）の質量速度である液体質量速度ＧLの関係について、計算結果と実験値との比較を示す図である。これらの図において、気体（Ｇ）の質量速度である気体質量速度ＧGがパラメータとしてある。これらの図において、計算結果は実線等で示され、実験結果はデータ点で示されている。図１２は、２つのピークを有する気孔率分布を１つの平均的気孔率分布で近似した場合の計算結果である。この場合、計算結果と実験値との間に大きな乖離があることが分かる。図１３は、図３で述べた手順で、複数の層にモデル化して計算した結果である。図１２に比較して、計算結果が実験値にかなりよく一致していることが分かる。
【図面の簡単な説明】
【００７４】
【図１】本発明に係る実施の形態における気液二相流れ評価シミュレーション装置の構成を示す図である。
【図２】気液二相流れの場合のダルシーの関係式を説明する図である。
【図３】本発明に係る実施の形態において、気液二相流れが供給されたときの多孔体の両端の圧力損失を算出する手順を示すフローチャートである。
【図４】本発明に係る実施の形態において、気孔率分布が複数のピークを有する例を説明する図である。
【図５】本発明に係る実施の形態において、気孔率分布が複数のピークを有する場合の多孔体の様子を説明する図である。
【図６】本発明に係る実施の形態において、モデル化を説明する図である。
【図７】本発明に係る実施の形態において、相対透過率算出の内部手順を示すフローチャートである。
【図８】本発明に係る実施の形態において、飽和度分布を求める様子を説明する図である。
【図９】本発明に係る実施の形態において、毛管圧分布を求める様子を説明する図である。
【図１０】本発明に係る実施の形態において、相対透過率分布を求める様子を説明する図である。
【図１１】本発明に係る実施の形態において、気体相対透過率と液体相対透過率の計算結果と、別途実験で求めた値との比較を示す図である。
【図１２】従来技術の方法で求めた多孔体の圧力損失と、実験値との比較を示す図である。
【図１３】本発明に係る実施の形態において、多孔体の圧力損失の計算結果と、実験値との比較を示す図である。
【符号の説明】
【００７５】
１０  気液二相流れ評価シミュレーション装置、１２  ＣＰＵ、１４  入力部、１６  出力部、１８  記憶装置、２０  モデル化モジュール、２２  相対透過率算出モジュール、２４  流入層ＰOUT1算出モジュール、２６  流出層ＰOUT算出モジュール、２８  圧力損失ΔＰ算出モジュール、３０  多孔体、３２  気液二相流れ。
【図１】
2008276357.tif 000003
【図２】
2008276357.tif 000004
【図３】
2008276357.tif 000005
【図４】
2008276357.tif 000006
【図５】
2008276357.tif 000007
【図６】
2008276357.tif 000008
【図７】
2008276357.tif 000009
【図８】
2008276357.tif 000010
【図９】
2008276357.tif 000011
【図１０】
2008276357.tif 000012
【図１１】
2008276357.tif 000013
【図１２】
2008276357.tif 000014
【図１３】
2008276357.tif 000015