JP4797344B2 - 冷媒熱物性値の計算方法、当該熱物性値の計算プログラム、当該計算プログラムを記録したコンピュータ読み取り可能な記録媒体及び冷媒熱物性値の計算装置 - Google Patents

Description

本発明は、冷媒熱物性値の計算に関するものであり、冷媒熱物性値を得るための高速で高精度な計算を実現し、また、解の安定性を実現した冷媒熱物性値の計算方法、当該熱物性値の計算プログラム、当該計算プログラムを記録したコンピュータ読み取り可能な記録媒体及びその計算プログラムを実装した計算装置に関するものである。

従来より、エアコンの設計を行うにあたって、構成要素である膨張弁、熱交換器、圧縮機の接続関係を指定し、さらには、熱交換器に用いられるフィンやチューブ等の形状をも指定して、冷凍サイクル全体の性能をシミュレートすることが出来るシミュレーション方法が提案されている。このような冷凍サイクルのシミュレーション方法を用いてエアコン等の設計を行うことにより、細部の変更等による全体への影響を実際の試作を行わずに知ることが出来るため、開発コストの削減等の効果を得ることが出来、また、開発に要する期間を短縮することが可能となる。このようなシミュレーション方法の一例としては、特許文献１及び特許文献２が既に提案されている。
特開平０９−２５７３１９号公報 特開平１１−２７０９１３号公報

このようなエアコンの設計においては、各部での冷媒の熱物性値を求めなければ全体の性能を予測することが出来ないため、冷凍サイクルのシミュレーションを行うにあたっては、冷媒熱物性を表現する関数が重要かつ不可欠である。冷凍サイクルのシミュレーションでは常に冷媒の熱物性値を計算すると同時に数値解析を行っているため、この関数は最終的に数千回も呼び出されて使用されることになる。よって、この冷媒熱物性の関数の１回の計算速度がシミュレーション全体の速度に大きな影響を与えることになる。

従来用いられてきた冷媒熱物性の関数として、Martin-Hou方程式を以下に示す。

上記（１）の式において、温度Ｔは定数であり、比体積νに対応する圧力ｐを計算する必要がある場合には、（１）式を変形して、次の陽方程式（２）を用いる。

圧力ｐに関しては（２）式の陽方程式を得ることが出来るため、比体積νに対応する圧力ｐを求める計算は非常に高速でかつ安定である。しかし、圧力ｐと温度Ｔに対応する比体積νを計算する場合、（１）式を変形しても、ν＝ｆ（ｐ）の形式となる陽方程式を得ることは不可能であるため、圧力ｐと温度Ｔに対応する比体積νを求めるには、陰方程式である（１）式に基づいて反復法で収束解を求めざるをえない。よって、計算速度は反復ループ数と反復回数の影響を受け、反復ループ数と反復回数が増えれば増えるほど計算速度は遅くなってしまう。また、反復演算の過程で発散が生じることで正しい解を得られない可能性もあり、このような場合には高速、高精度、高安定というシミュレーションに求められる要請に答えることができない。

冷媒熱物性のうち、温度Ｔ、比体積ν及び圧力ｐの関係について（１）式を用いて説明したが、この他の冷媒熱物性として、例えば、温度Ｔ、比体積ν、圧力ｐ、エンタルピｈ、エントロピｓ、密度ρ、定圧比熱Ｃ_ｐ、動粘性μ、熱伝導率λ、表面張力σ、乾度χ等のパラメータについても相互の関係を表した関数が存在し、それらの関数から陽方程式を導出できない限り、比体積νと圧力ｐの関係と同様に、陰方程式に基づいて反復法で収束解を求めるしかなく、この場合においても、反復ループ数と反復回数が増えれば増えるほど計算速度は遅くなるし、また、反復演算の過程で発散が生じるおそれもある。

このように、従来の冷媒熱物性の計算方法は、冷媒熱物性の算出方法として、複雑な関数において反復計算による収束解を求める方法で冷媒の状態量から冷媒の熱物性値を求めている。この方法は、冷媒を使用する可能性のあるあらゆる条件範囲を広くカバーしており、また、高精度で冷媒熱物性を計算できるという長所を有しているが、上述の通り複雑な関数による演算は非常に長い演算時間を要し、また、収束しない可能性もある不安定な計算方法であった。エアコンを設計するにあたって、高精度で冷媒熱物性を計算できることは大きなメリットであるが、シミュレーションの度に膨大な計算時間を必要とする従来方法は作業効率が悪く、シミュレータとしての実用性は著しく低いものとなってしまっていた。

本発明は、上記問題点に鑑みなされたものであり、冷媒熱物性値を得るための高速で高精度な計算を実現し、また、解の安定性を実現することで、シミュレータとしての実用性を向上させた冷媒熱物性値の計算方法、当該熱物性値の計算プログラム、当該計算プログラムを記録したコンピュータ読み取り可能な記録媒体及びその計算プログラムを実装した計算装置を提供することを目的とするものである。

本発明の請求項１は、冷媒熱物性値のソースデータを得るためのソースデータ演算処理手段と、特定の高次多項式を前記冷媒熱物性値のソースデータのうち、求めたい冷媒熱物性値によって回帰処理することで前記高次多項式の係数を得る回帰処理手段と、前記求めた係数を代入した高次多項式を解くことで前記高次多項式の根を得て、前記根の中から有効な根を選択する解式及び根選択処理手段とによって、求めたい冷媒熱物性値の関係を表した陽方程式を得て、この陽方程式を用いて冷媒熱物性値の計算を行うことを特徴とする冷媒熱物性値の計算方法である。

本発明の請求項２は、請求項１に加えて、ソースデータ演算処理手段において求めるソースデータは、冷媒熱物性値の関係を表した陰方程式を反復法によって解いて求めたものを含むことを特徴とする冷媒熱物性値の計算方法である。

本発明の請求項３は、請求項１又は２に加えて、ソースデータ演算処理手段では、指定した冷媒の温度範囲についてソースデータを求め、回帰処理手段では、前記冷媒の温度範囲を指定して求めたソースデータのうち、求めたい冷媒熱物性値について高次多項式に当てはめて最小二乗法により回帰処理し、この回帰処理が収束した場合に前記高次多項式の係数を得るようにしたことを特徴とする冷媒熱物性値の計算方法である。

本発明の請求項４は、請求項３に加えて、回帰処理手段において行った回帰処理が収束しなかった場合には、冷媒の温度範囲を調整し、その調整後の温度範囲のソースデータを用いて回帰処理を行うようにし、温度範囲を調整しながら回帰処理を設定した回数行っても収束しない場合には、回帰処理に用いたソースデータに対して変換式を適用して再度回帰処理を行うようにしたことを特徴とする冷媒熱物性値の計算方法である。

本発明の請求項５は、請求項１に加えて、解式及び根選択処理手段では、求めた前記高次多項式の根をそれぞれソースデータと比較して、最も誤差の小さい根を有効な根とすることを特徴とする冷媒熱物性値の計算方法である。

本発明の請求項６は、請求項１乃至５記載のソースデータ演算処理手段、回帰処理手段、及び、解式及び根選択処理手段とをコンピュータに実行させ、これによって求めたい冷媒熱物性値の関係を表した陽方程式を得て、この陽方程式を用いて冷媒熱物性値の計算を行うことを特徴とする冷媒熱物性値の計算プログラムである。

本発明の請求項７は、請求項１乃至５記載のソースデータ演算処理手段、回帰処理手段、及び、解式及び根選択処理手段とをコンピュータに実行させ、これによって求めたい冷媒熱物性値の関係を表した陽方程式を得て、この陽方程式を用いて冷媒熱物性値の計算を行うことを特徴とする冷媒熱物性値の計算プログラムを記憶した記憶媒体である。

本発明の請求項８は、冷媒熱物性値のソースデータを得るためのソースデータ演算処理部と、特定の高次多項式を前記冷媒熱物性値のソースデータのうち、求めたい冷媒熱物性値によって回帰処理することで前記高次多項式の係数を得る回帰処理部と、前記求めた係数を代入した高次多項式を解くことで前記高次多項式の根を得て、前記根の中から有効な根を選択する解式及び根選択処理部とを具備し、これらの処理部を経ることで求めたい冷媒熱物性値の関係を表した陽方程式を得てこれをメモリに記憶し、この陽方程式を用いて冷媒熱物性値の計算を行うことを特徴とする冷媒熱物性値の計算プログラムを実装した冷媒熱物性値の計算装置である。

本発明の請求項９は、請求項８に加えて、ソースデータ演算処理部において求めるソースデータは、冷媒熱物性値の関係を表した陰方程式を反復法によって解いて求めたものを含むことを特徴とする冷媒熱物性値の計算プログラムを実装した冷媒熱物性値の計算装置である。

本発明の請求項１０は、請求項８又は９に加えて、ソースデータ演算処理部では、指定した冷媒の温度範囲についてソースデータを求め、回帰処理部では、前記冷媒の温度範囲を指定して求めたソースデータのうち、求めたい冷媒熱物性値について高次多項式に当てはめて最小二乗法により回帰処理し、この回帰処理が収束した場合に前記高次多項式の係数を得るようにしたことを特徴とする冷媒熱物性値の計算プログラムを実装した冷媒熱物性値の計算装置である。

本発明の請求項１１は、請求項１０に加えて、回帰処理部において行った回帰処理が収束しなかった場合には、冷媒の温度範囲を調整し、その調整後の温度範囲のソースデータを用いて回帰処理を行うようにし、温度範囲を調整しながら回帰処理を設定した回数行っても収束しない場合には、回帰処理に用いたソースデータに対して変換式を適用して再度回帰処理を行うようにしたことを特徴とする冷媒熱物性値の計算プログラムを実装した冷媒熱物性値の計算装置である。

本発明の請求項１２は、請求項８に加えて、解式及び根選択処理部では、求めた前記高次多項式の根をそれぞれソースデータと比較して、最も誤差の小さい根を有効な根とすることを特徴とする冷媒熱物性値の計算プログラムを実装した冷媒熱物性値の計算装置である。

上記のような冷媒熱物性値の計算方法、計算プログラム、計算プログラムを記憶した記憶媒体、又は、計算プログラムを実装した冷媒熱物性値の計算装置とすることで、従来計算に長時間を要し、さらに解が発散する可能性のあった反復法による演算に替えて、冷媒熱物性値を得るための高速で高精度な計算を実現し、また、解の収束の安定性を実現することが可能な陽方程式を得ることが出来たので、これらの陽方程式を用いることで実用性の高いシミュレーションを実現できる。

本発明による冷媒熱物性値の計算は、陰方程式を解くことでしか求めることができなかった冷媒熱物性のソースデータに対して回帰処理を行うことで近似陽方程式を求め、この近似陽方程式を用いて冷媒熱物性値の計算を行うことで、高速で高精度な計算を実現し、また、解の安定性を実現しているものである。この本発明の冷媒熱物性値の計算を実現するための全体の流れを図１のフローチャート図を用いて説明する。この図１のフローチャートに示すように、本発明は３つの主要な処理系として、ソースデータ演算処理系（S102〜S103）、回帰処理系（S104〜S108）、解式及び根選択処理系（S109〜S111）を有している。以下、３つの処理系の概要を説明する。

［ソースデータ演算処理系］
先ず、（S101）においてｋ＝１とした後に、（S102）において、ソースデータを演算する必要のあるデータ範囲を特定するために、最低温度Ｔ_min、最高温度Ｔ_max及びこれらの温度間で計算するデータ点の個数Ｎを指定する。ここで指定する温度範囲は、冷媒の種類によって異なるものであり、各冷媒の具体的な温度範囲については後述する。温度範囲を指定したら、（S103）においてソースデータを演算によって求める。ここで求めるソースデータは、従来の反復法を用いた計算方法によって求められるものであり、例えば、温度Ｔ、比体積ν、圧力ｐ、エンタルピｈ、エントロピｓ、密度ρ、定圧比熱Ｃ_p、動粘性μ、熱伝導率λ、表面張力σ、乾度χ等の冷媒熱物性値を全てソースデータとして演算によって求めてメモリ等に記憶する。

［回帰処理系］
次に、前記（S103）において演算によって求めたソースデータのうち、任意のパラメータについての冷媒熱物性値を（S104）で読込み、この読込んだソースデータを（S105）において回帰式に入力するのに最適な形式となるように変換処理を行う。そして、（S106）において変換処理後のソースデータを回帰式に入力して、（S107）において回帰処理を行う。ここで、（S106）において使用する回帰式は、多次多項式を変形したものであり、未知変数部分にそれぞれ変換処理後のソースデータを入力して回帰処理を行うことにより、前記多次多項式の係数が求まる。この多次多項式は冷媒の種類（純冷媒、混合冷媒など）やその状態（サブクール領域、飽和領域又はスーパーヒート領域）によって異なり、適宜選択して使用する。このようにして回帰処理を行った後に、その回帰処理が収束したか否かを（S108）において判断する。回帰処理が収束した場合には次の（S109）に進むが、回帰処理が収束しない場合には、（S116）において指定した温度範囲及びその間のデータ数を調整し直して、再度（S103）〜（S107）によって回帰処理を行って、（S108）で収束判断を行う。回帰処理を３回行っても収束しない場合（（S113）においてnoの場合）には、（S114）で変換式を適用することで、ソースデータが回帰式に適用可能な形式になるように（S105）で変換処理を行う。

［解式及び根選択処理系］
回帰処理が収束した場合には、前記多次多項式の係数が求まるのでこれを（S109）でメモリ等に記憶し、これらの係数を用いて残りの係数を全て求めて（S110）でメモリ等に記憶し、全ての係数が求まったら、これらを（S110）において前記多次多項式に代入してこれを解く。根の数は多項式の次数で決まり、これらの根の中から最適な根を選択する。

このように、ソースデータ演算処理系、回帰処理系、解式及び根選択処理系の３つの処理系によって、任意に選択したパラメータについての関係を近似式としての陽方程式によって表現できるようになる。この陽方程式によって冷媒熱物性値を求める場合には、反復計算を必要としないため、高い計算安定性を有するとともに、高精度でかつ高速に計算を行うことが可能となる。また、この方法で開発した陽方程式のうちのいくつかは、冷媒熱物性の計算における可逆性を備えているものである。
上記の処理系の説明では概要のみを記載したが、実際には冷媒の種類（純冷媒、混合冷媒）や冷媒の状態（サブクール領域、飽和領域又はスーパーヒート領域）によって高速計算の実現方法はそれぞれ異なる。

具体的な条件時における高速計算の実現方法を説明する前に、本発明の実施例において使用する冷媒の種類とその温度範囲を説明する。純冷媒としては、Ｒ２２（飽和温度範囲は−６０℃〜＋８０℃）、Ｒ１３４ａ（飽和温度範囲は−６０℃〜＋８０℃）、Ｒ１２５（飽和温度範囲は−６０℃〜＋６０℃）、Ｒ３２（飽和温度範囲は−６０℃〜＋７０℃）などがあり、これらは全てスーパーヒート（過熱度）を＋６５℃までとする。混合冷媒としては、Ｒ４１０Ａ（飽和液温度範囲は−６０℃〜＋６０℃、飽和蒸気温度範囲は−６０℃〜＋６０℃）、Ｒ４０７Ｃ（飽和液温度範囲は−６０℃〜＋８０℃、飽和蒸気温度範囲は−６０℃〜＋８０℃）などがあり、これらもスーパーヒート（過熱度）を＋６５℃までとする。これらの冷媒について、以下、それぞれの具体的な条件時における高速計算の実現方法を説明する。

［純冷媒が飽和領域にある場合］
純冷媒としては、Ｒ２２、Ｒ１３４ａ、Ｒ１２５、Ｒ３２などが存在するが、これらの純冷媒が飽和領域にある場合の冷媒熱物性を求めるための陽方程式を導出するのに、次の３次方程式を用いた。

ここで、ｕ、ｖは、２つの異なる温度特性を表した変数であり、ａ_１〜ａ_９は、係数である。この式（３）において一方の温度特性ｖの値が利用可能である場合には、式（３）は以下の式（４）に変換することが出来る。

この式（４）は３次方程式であるので、一般に、これは３つの異なる根ｕ_１、ｕ_２及びｕ_３を持つことになり、それは次のように表現される。

これらの式（５）、式（６）、式（７）のうち、ただ１つの式が必要とする正しい関数であり、これは３つの方程式を計算してその結果をソースデータと比較することで正しい関数を選択することができる。このようにして選択された関数は、ｕとｖの正しい関係を表した反復計算の必要がない陽方程式となるため、高速演算と安定性を保証できるものとなる。

また、上記式（３）において他方の温度特性ｕの値が利用可能である場合には、式（３）は以下の式（８）に変換することが出来る。

この式（８）は３次方程式であるので、一般に、これは３つの異なる根ｖ_１、ｖ_２及びｖ_３を持つことになり、それは次のように表現される。

これらの式（９）、式（１０）、式（１１）のうち、ただ１つの式が必要とする正しい関数であり、これは３つの方程式を計算してその結果をソースデータと比較することで正しい関数を選択することができる。このようにして選択された関数は、ｖとｕの正しい関係を表した反復計算の必要がない陽方程式となるため、高速演算と安定性を保証できるものとなる。
上記式（４）と式（８）は同じ方程式（３）を変形したものであり、式（４）の有効な根が式（６）で、式（８）の有効な根が式（９）であるとすると、これらの式から得られたｕとｖは一対一の対応をしており、それぞれ一方の式から他方の根を逆算できる可逆性を備えたものとなる。

よって、冷媒熱物性値のうち選択した２つのパラメータについてのソースデータを用いて式（３）を回帰処理して係数ａ_１〜ａ_９を求めることで、選択した２つのパラメータについての陽方程式（例えば、式（６）や式（９））を導出することが出来る。ここで、式（３）の回帰処理を簡単にするために、回帰処理の前に式（３）に適切な処理をする必要がある。あるデータ範囲において、点（ｘ_０、ｙ_０）、点（ｘ_１、ｙ_１）及び点（ｘ、ｙ）は、それぞれ最低点、最高点及びランダムな点であると仮定して、ｕ＝ｙ−ｙ_０、ｖ＝ｘ−ｘ_０と置くと、式（３）は次のようになる。

２つの端点（ｘ_０、ｙ_０）及び（ｘ_１、ｙ_１）について、式（１２）は次のようになる。

この式（１３）と式（１４）とから、次の式（１５）及び式（１６）が得られる。

式（１６）はデータ範囲の全ての点に適用可能なので、式（１５）と式（１６）を式（３）に代入すると、次の式を得ることができる。

この式（１７）を用いて最小二乗法により選択した前記データ範囲の全データ点について回帰処理を行うことにより、７個の未知係数ａ_１〜ａ_７を得ることができ、この係数ａ_１〜ａ_７を式（１６）に代入することで未知係数ａ_８を得ることができる。よって、この式（１７）は純冷媒の飽和領域における冷媒熱物性を求めるための陽方程式を導出するための回帰式となる。

上記回帰式（１７）を用いた回帰処理によって、純冷媒が飽和領域にある場合の冷媒熱物性値を求めるための陽方程式を導出する流れを図１のフローチャートを用いて詳しく説明する。
先ず、（S101）において回帰処理の回数ｋをｋ＝１とした後に、（S102）において、ソースデータを演算する必要のあるデータ範囲を特定するために、最低温度Ｔ_min、最高温度Ｔ_max及びこれらの温度間で計算するデータ点の個数Ｎを指定する。例えば、冷媒がＲ２２である場合には、有効なデータ範囲は−６０℃〜＋８０℃であり、この間のデータ点の間隔を１℃とするとデータ点の個数Ｎは１４１個となる。この他の純冷媒の有効なデータ範囲としては、Ｒ１３４ａの場合は−６０℃〜＋８０℃、Ｒ１２５の場合には−６０℃〜＋６０℃、Ｒ３２の場合には−６０℃〜＋７０℃となっている。

温度範囲を指定したら、（S103）においてソースデータを演算によって求める。ここでのソースデータ演算処理の詳しいフローチャートを図２に示す。ここで、ｉをデータ点の番号（ｉ＝１、２、・・・、１４１）とする。この図２において、（S201）でｉ＝１に設定した後に、（S202）において、Ｔ_s［ｉ］＝ｉ＋Ｔ_smin−１の式によってｉ番目のデータ点の温度を指定する。ここで、Ｔ_sminは飽和領域における最低温度を表している。（S202）で最初の温度を指定したら、（S203）において温度Ｔ_s［ｉ］における冷媒熱物性の様々なパラメータを演算によって求める。ここで求められるソースデータは、従来の反復法を用いた計算方法によって求められるもので、例えば、液相の圧力ｐ_l［ｉ］、気相の圧力ｐ_v［ｉ］、液相のエンタルピｈ_l［ｉ］、気相のエンタルピｈ_v［ｉ］等があり、この他の冷媒熱物性値についても全てソースデータとして演算によって求める。この最初の（S202）と（S203）のフローによって最低温度Ｔ_sminにおける全ての冷媒熱物性値を演算によって求めて、次の（S204）においてｉ＜Ｎである場合には（S205）においてｉの値に＋１をして、以下、Ｎ個分のデータ点について冷媒熱物性値を求めるまで、（S202）〜（S205）の処理を続ける。ｉ＝Ｎになったら、（S206）においてＮ個分のデータ点のそれぞれにおいて演算で求めた全ての冷媒熱物性値のソースデータをメモリ等に記憶して、図1の（S104）に進む。

（S104）では、前記（S103）において演算によって求めたソースデータのうち、任意の２つのパラメータについての冷媒熱物性値を読込み、この読込んだソースデータを（S105）において回帰式に入力するのに最適な形式となるように変換処理を行う。ここで、（S104）において読込むソースデータを圧力ｐ_l［ｉ］と温度Ｔ_s［ｉ］とした場合の（S105）でのソースデータの処理の詳細を図３のフローチャートを用いて説明する。この図３において、先ず（S301）においてｋ＝３か否かを確認し、ｋ＝３でない場合には（S302）に進む。この（S302）では、回帰式（１７）を回帰処理する際に用いる値として、ｕ［１］＝ｐ_l［Ｎ］−ｐ_l［１］、ｖ［１］＝Ｔ_s［Ｎ］−Ｔ_s［１］の２つを準備する。次に（S303）においてｉ＝２とした後に、（S304）〜（S306）のループによってｉ＝２〜ｉ＝Ｎとなるまでの全データ点について、ｕ［ｉ］＝ｐ_l［ｉ］−ｐ_l［１］、ｖ［ｉ］＝Ｔ_s［ｉ］−Ｔ_s［１］の２つを準備する。ｉ＝Ｎまで２つの変数ｕ［ｉ］とｖ［ｉ］を準備したら、（S305）で判断して図１の（S106）へ進む。

ここで、前記（S301）においてｋ＝３であった場合、すなわち回帰処理を３回行ったが収束しなかった場合には、図１の（S114）において回帰処理が収束するように変換式を適用する。この図３の場合には、（S301）にてｋ＝３であった場合には（S307）に進み、この（S307）においてｉ＝1とした後に、前記（S104）において読込んだ圧力ｐ_l［ｉ］とＴ_s［ｉ］のうち、圧力ｐ_l［ｉ］についてのみ（S308）においてｙ´＝ｙ／１００という変換式を適用することにより１００分の１にする。これをｉ＝１〜Ｎまで、すなわち読込んだ全ての圧力ｐ_l［ｉ］について、（S308）〜（S310）のループによって変換式を適用して、変換終了後に、前記（S302）からの処理に移行する。ここでの変換式の適用は、回帰計算の解を容易に求めるため、ソースデータに対する一組の変数間の関係を調整するために行う。これらの変換式は、実験的、経験的に決定される。

前記（S105）でのソースデータの処理が終了したら、その処理後のｕ［ｉ］とｖ［ｉ］とを（S106）において回帰式（１７）に入力し、（S107）においてその回帰式（１７）を回帰処理する。この（S107）において行う回帰処理は、一般的な数学的手法として知られている最小二乗法により行われるものであり、具体的な回帰処理の内容については説明を省略する。この回帰処理によりソースデータに対して最適な曲線を適用することができ、その曲線の係数としてａ_１〜ａ_７を得ることができる。ここでの回帰処理が収束したか否かを（S108）において判断するが、収束しなかった場合には、（S115）でｋの値をカウントアップしつつ（S116）で最低温度Ｔ_min、最高温度Ｔ_max及びデータ点の個数Ｎを調整し直して、再度（S103）〜（S107）の処理を行う。温度範囲の指定があまりにも広すぎると収束しない可能性があるため、最適な温度範囲に調整することで回帰処理を収束させることが可能となる。ただし、温度範囲を調整しても回帰処理が収束しない場合もあり、そのような場合には、前記（S105）でのソースデータの処理の部分で説明した通り、３回の回帰処理で収束しない場合に限り（S114）で変換式を適用することとなる。

（S108）で回帰処理が収束したら、（S109）で前記回帰処理によって得られた係数ａ_１〜ａ_７をメモリ等に記憶し、（S110）で係数ａ_１〜ａ_７を式（１６）に代入することで係数ａ_８を得てこれをメモリ等に記憶する。このようにして得られた係数ａ_１〜ａ_８を用いて（S111）で式（４）及び式（８）を解いて、それらの根の中から最適な根を選択する。この（S111）における処理の詳細を図４に基づいて説明する。図４において、（S401）で式（４）に係数ａ_１〜ａ_８を代入し、（S402）で式（４）を解く。式（４）を解くと複数の根が求められるので、（S403）でこれらの根の中から有効な根を選択する。次に、（S404）で式（８）に係数ａ_１〜ａ_８を代入し、（S405）で式（８）を解く。式（８）を解くと複数の根が求められるので、（S406）でこれらの根の中から有効な根を選択する。最後に、（S407）で式（４）の有効な根と式（８）の有効な根とをそれぞれ出力する。この式（４）の有効な根と式（８）の有効な根とが、目的とする冷媒熱物性値を求めるための陽方程式となる。

上記（S402）で式（４）を解く具体的な流れを図５のフローチャートに基づいて説明する。図５において、（S501）でｉ＝１とした後、（S502）で回帰処理が収束するのに要した回数としてｋ＝３であるか否かを確認する。３回以下で回帰処理が収束している場合にはそのまま（S504）へ進むが、回帰処理が３回以下で収束していない場合には、（S114）を経て（S308）で変換式を適用しているので、それに合わせて（S503）で最低温度における圧力値を１００分の１に変換する。次に、（S504）でｘ＝Ｔ_s［ｉ］−Ｔ_s［１］と定義し、また、（S505）において、前記ｘの値と（S110）で得られた係数ａ_１〜ａ_８を図９に示す式（Ａ−６）に適用することでａ、ｂ、ｃ及びｄの値を求める。この求めたａ、ｂ、ｃ及びｄの値を用いて（S506）において、図８の式（Ａ−１）〜式（Ａ−３）、及び、図９の式（Ａ−４）と式（Ａ−５）の計５つの式によって、ｙ₁［ｉ］〜ｙ₅［ｉ］を得る。ここで、式（Ａ−１）〜式（Ａ−３）は、式（４）の判別式が０より小さい場合の３つ実根であり、式（Ａ−４）と式（Ａ−５）は、式（４）の判別式が０より大きい場合の１つの実根と２つの虚根の対である。次に、（S507）において、ｙ₁［ｉ］〜ｙ₅［ｉ］にｐ_l［１］をそれぞれ加えることでｐ₁［ｉ］〜ｐ₅［ｉ］を得る。ここで、変換式を適用していた場合には、（S509）においてｐ₁［ｉ］〜ｐ₅［ｉ］をそれぞれ１００倍する。このような処理を（S510）においてｉがデータ点の個数Ｎとなるのを確認するまで、（S511）でｉを１ずつ増やしながらｐ₁［ｉ］〜ｐ₅［ｉ］をその都度求める。よって、この図５のフローチャートによる処理を終了した段階で、（ｐ₁［ｉ］〜ｐ₅［ｉ］）×Ｎ個分のデータが得られている。なお、純冷媒Ｒ２２における圧力ｐ_l［ｉ］と温度Ｔ_s［ｉ］の関係については、３つの実根を持つので、以下、式（Ａ−１）〜式（Ａ−３）で求めた（ｐ₁［ｉ］〜ｐ₃［ｉ］）×Ｎ個分のデータに基づいて説明する。

上記（S402）で式（４）を解くことで得た（ｐ₁［ｉ］〜ｐ₃［ｉ］）×Ｎ個分のデータを用いて、（S403）で有効な根を選択する処理について、図６のフローチャートを用いて説明する。有効な根を選択する手法は、反復法で求めたソースデータとしての圧力ｐ_l［ｉ］に対して、式（Ａ−１）〜式（Ａ−３）の３つの式でそれぞれ求めたｐ₁［ｉ］〜ｐ₃［ｉ］のどれが最も誤差が少ないかを判断の材料とする。先ず、（S601）においてｉ＝１とした後に、（S602）において、δ₁［ｉ］＝｜ｐ₁［ｉ］−ｐ_l［ｉ］｜／ｐ_l［ｉ］というようにしてδ₁［ｉ］を求め、同様に、δ₂［ｉ］、δ₃［ｉ］についても求める。この、δ₁［ｉ］〜δ₃［ｉ］は、式（Ａ−１）〜式（Ａ−３）で求めたデータの各点でのソースデータとの誤差である。これを（S602）〜（S604）の処理をｉの値を１ずつ増やしながら繰り返すことで、Ｎ個のデータ点全てについて各点での誤差を求める。次に、（S605）において、各点の誤差値であったδ₁［ｉ］〜δ₃［ｉ］について１〜Ｎまでを合計したものとして、δ₁、δ₂、δ₃を求める。（S606）において、このδ₁、δ₂、δ₃のうち最も値の小さいものをδとして選択する。例えば、ここでδ₂が最小であったとすると、次の（S607）において、δ_mean＝δ₂／Nによって誤差の平均を求め、δ_max＝ｍａｘ（δ₂［ｉ］）によって最大の誤差となったデータ点を求める。そして、（S608）において、δ_mean＜0.002、かつ、δ_max＜0.2という条件を満たした場合には、式（Ａ−２）で求めた根を有効な根として選択する。

このように図６のフローチャートによって、（S403）で式（４）についての有効な根を選択することができる。また、（S404）〜（S406）によって、式（８）を解いてその解の中から有効な根を選択するが、この場合の具体的なフローについても、図５及び図６で示した式（４）に関するものと略同様であり、図５の（S503）において用いた変換式を、図７の（S701）に置換え、図５及び図６のフローチャート中の「ｐ_ｌ」と「Ｔ_ｓ」とを置換え、また、図５において使用した式（Ａ−６）を図９に示す式（Ａ−７）に置換ることで、式（８）を解いてその解の中から有効な根を選択することができる。そして、図４のフローの最後として（S407）において、選択した式（４）の有効な根、及び、式（８）の有効な根をそれぞれ出力する。

このような流れにより、図１の（S111）において、式（４）の有効な根、及び、式（８）の有効な根を得ることができ、これらは、ｐ_ｌ＝ｆ（Ｔ）、及び、Ｔ_ｓ＝ｆ（ｐ_ｌ）というように表現された圧力と温度の関係を表した陽方程式である。この陽方程式を用いることにより、未知の飽和液圧力ｐ_ｌ又は未知の飽和液温度Ｔ_ｓを得るための、高速で高精度かつ安定な計算を実現することができる。このように、図１のフローチャートによって、他のパラメータについても冷媒熱物性値を求めるための陽方程式を導出することができ、その導出した陽方程式を使った計算は高速で高精度かつ高安定なものとなる。

［純冷媒がスーパーヒート領域にある場合］
前記実施例１では、純冷媒が飽和領域にある場合の冷媒熱物性を求めるための陽方程式を導出する方法について説明したが、本実施例２では、純冷媒がスーパーヒート領域にある場合の冷媒熱物性を求めるための陽方程式を導出する方法について説明する。この純冷媒がスーパーヒート領域にある場合についての冷媒熱物性を求めるための陽方程式を導出するのに、次の３次方程式を用いた。

ここで、ｙはスーパーヒート温度特性、ｙ_ｓはｙに対応した飽和蒸発温度特性、Ｔはスーパーヒート温度、Ｔ_ｓはＴに対応した飽和蒸発温度であるとともにスーパーヒート圧力ｐの関数であり、ａ_１〜ａ_１１は係数である。この式（１８）は、それぞれｙが未知か、Ｔが未知かの何れかによって、式（１９）又は式（２０）に変換することが出来る。

この式（１９）及び式（２０）は、共に３次方程式であるので、一般に、これらは３つの異なる根を持つことになる。それらの根のうち、式（１９）と式（２０）はそれぞれ１つの有効な根を持ち、それらは、例えばｙ＝ｆ（ｐ，Ｔ）のように、陽方程式として表現することが可能となる。このように表現された陽方程式に反復計算は必要ないため、高速で安定な計算を行うことができる。

よって、上記式（１９）及び式（２０）を回帰処理することにより、純冷媒がスーパーヒート領域にある場合の冷媒熱物性を求めるための陽方程式を導出することができる。この場合の流れは、いくつかのステップを調整する以外は、図１のフローチャートと基本的に同様である。以下、調整したステップの部分について説明する。

前記実施例１では、図１の（S102）において、ソースデータを演算する必要のあるデータ範囲を特定するために、最低温度Ｔ_min、最高温度Ｔ_max及びこれらの温度間で計算するデータ点の個数Ｎを指定し、回帰処理が収束しない場合には（S116）において前記最低温度Ｔ_min、最高温度Ｔ_max及びデータ点の個数Ｎを調整するようになっていたが、純冷媒がスーパーヒート領域にある場合には、この図１の（S102）を図１０（ａ）に示す（S1001）に置換えて、同様に、図１の（S116）を図１０（ｂ）に示す（S1002）に置換える。ここで、Ｔ_sminは最低飽和蒸発温度、Ｔ_smaxは最高飽和蒸発温度、Ｎはこれらの温度間で計算するデータ点の個数、Ｎ_spはスーパーヒート（過熱度）であり、これらを（S1002）において指定してから、次の（S103）に進む。同様に、回帰処理が収束しない場合には（S1002）において最低飽和蒸発温度Ｔ_smin、最高飽和蒸発温度Ｔ_smax、これらの温度間で計算するデータ点の個数Ｎ、及び、スーパーヒート（過熱度）Ｎ_spを調整する。例えば、純冷媒がＲ２２の場合には、最低飽和蒸発温度Ｔ_smin＝−６０℃、最高飽和蒸発温度Ｔ_smax＝＋８０℃、これらの温度間で計算するデータ点の個数Ｎ＝１４１個、及び、スーパーヒート（過熱度）Ｎ_sp＝＋６５℃となる。

前記実施例１では、図１の（S103）においてソースデータを演算によって求める詳しいフローチャートとして図２を用いたが、純冷媒がスーパーヒート領域にある場合には、この図２の替わりに図１１を用いる。図１１に示すソースデータの演算の流れにおいて図２と異なる点は、図２においては、ある温度Ｔ［ｉ］におけるソースデータを求めてこれをＮ点分求めるものであったが、純冷媒がスーパーヒート領域にある場合には、図１１の（S1104）〜（S1108）のループにおいて、ある飽和蒸発温度Ｔ_s［ｉ］においてスーパーヒート（過熱度）温度を１度ずつ上げながらＮ_spに到達するまで各スーパーヒート温度毎にソースデータを演算し、この処理を（S1110）のループによって飽和蒸発温度Ｔ_s［ｉ］を１度ずつ上げながらNに到達するまで繰り返す。よって、純冷媒Ｒ２２がスーパーヒート領域にある場合には、（ソースデータのパラメータ数×６５×１４１）個のソースデータを演算によって求めることとなる。

このようにしてソースデータを演算によって求めたら、図１のフローチャートの（S104）〜（S111）の流れと同様に、求めたソースデータを用いて式（１８）を回帰処理をし、回帰処理によって求まった係数を用いて式（１８）を解いて有効な根を選択する。ここでの回帰処理、並びに、解式及び根の選択処理は、若干の変換式の違いがある以外は、前記実施例１での最小二乗法を用いた回帰処理と略同様であり、また、根の選択方法についても前記実施例１と略同様であるので、説明は省略する。

［飽和領域又はスーパーヒート領域で、いくつかの特別冷媒温度特性を求める場合］
計算過程では、未知変数、また稀ではあるが既知変数としてのみ用いるいくつかの特別冷媒温度特性として、例えば、定圧比熱Ｃp、動粘性μ、熱伝導率λ、表面張力σなどがある。例えば、既知の飽和温度Ｔ_sによって未知の飽和熱伝導率λを計算することが多いし、稀には、既知の飽和熱伝導率λによって未知の飽和温度Ｔ_sを計算することもある。このような場合、この種の温度特性では、一組のλとＴ_sの一対一関数を得るために３次方程式を解いて陽関数を得る必要はなく、高次陽多項式を適用して、高速、高精度かつ高安定に計算できる冷媒熱物性の陽方程式を得ることができる。ここで、飽和領域（サブクール領域を含む）又はスーパーヒート領域におけるいくつかの特別冷媒温度特性を求めるために、以下の高次陽多項式を用いる。

この式（２１）に対して求めたい特別冷媒温度特性のソースデータ、例えば、λとＴ_sのソースデータを適用して回帰処理を行うが、その回帰処理は、若干の変換式の違いがある以外は、冷媒が飽和状態かスーパーヒートの常態かによって、前記実施例１で述べた最小二乗法によるものと略同様又は前記実施例２で述べた最小二乗法によるものと略同様であるので、説明は省略する。

［純冷媒の二相冷媒温度特性を求める場合］
二相純冷媒の場合、圧力ｐ_tp（或いは温度Ｔ_tp）を与えたとき、任意の気相温度特性ｚ_vと液相温度特性ｚ_lは、それぞれ対応する飽和蒸気温度特性ｚ_vsと飽和液温度特性ｚ_lsに等しくなる。乾度χが０〜１の間の領域においては、対応する飽和蒸気温度特性ｚ_vs並びに飽和液温度特性ｚ_lsは変化しない。同時に、気相冷媒の温度特性は対応する飽和蒸気冷媒に等しく、液相冷媒の温度特性は対応する飽和液冷媒に等しくなる。その結果、二相純冷媒は、圧力（或いは温度）及びその他の異なる温度常態パラメータ、例えば、乾度χ、エントロピｓ、エンタルピｈなどを与えると決めることができる。よって、例えば、圧力ｐ_tp及び乾度χを与えると、任意の対応する二相温度特性ｚ_tpは、次式で計算できる。
ｚ_tp＝（１−χ）ｚ_ls+χｚ_vs （22）

この式（２２）において、飽和蒸気温度特性ｚ_vsと飽和液温度特性ｚ_lsとが明らかになれば二相温度特性ｚ_tpは求まるので、飽和蒸気温度特性ｚ_vsと飽和液温度特性ｚ_lsとをそれぞれ前記実施例における方法で求めることにより、二相温度特性ｚ_tpが得られる。

［混合冷媒の二相冷媒温度特性を求める場合］
混合冷媒（非共沸冷媒）として、例えば、Ｒ４１０Ａ、Ｒ４０７Ｃなどが存在する。混合冷媒、特に非共沸冷媒の場合、二相圧力ｐ_tpが与えられたとき、二相温度特性ｚ_tpは、対応する沸点温度特性ｚ_bにも等しくないし、対応する露点温度特性ｚ_dにも等しくない。同時に純冷媒とも異なり、ある圧力に対応する沸点温度Ｔ_bと対応する露点温度Ｔ_dも等しくない。従って、混合冷媒の場合、気相状態は対応する露点温度状態に等しくなく、また、液相状態は対応する沸点温度状態に等しくもない。また、純冷媒とは異なり、混合冷媒の温度特性は対応する飽和（沸点或いは露点）温度特性を直接用いて計算することができない。ここで、液相温度特性（沸点温度特性ではない）をｚ_lとし、気相温度特性（露点温度特性ではない）をｚ_vとしたとき、対応する任意の温度特性ｚ_tpは、次式で計算できる。
ｚ_tp＝（１−χ）ｚ_l+χｚ_v （23）

この式（２３）をｚ_lとｚ_vとによって回帰することで、陽方程式を得られるように見えるが、液相温度特性ｚ_lと気相温度特性ｚ_vは、それぞれ液相冷媒と気相冷媒の要素に準拠しており、これらの要素は乾度χによって変化するため、各乾度χにおける液相温度特性ｚ_lと気相温度特性ｚ_vを計算するのに全ての高速演算式を回帰するのは不可能である。よって、混合冷媒の温度特性を求めるために次の新たな式を導入した。

一般に、混合冷媒の圧力とエンタルピ或いは圧力と乾度を用いることができる場合、そのときの物理状態は決まる。よって、混合冷媒の計算を行うためには、圧力、温度、エンタルピ及びエントロピのみが必要である。この条件に基づいて、混合冷媒の温度特性を求めるために次の４つの式を用いる。１つは、圧力と乾度が既知である場合にエンタルピを計算するための式で、他の３つは、圧力とエンタルピが既知である場合に、乾度、エントロピ及び温度を計算するための式である。

以下の式において、エンタルピｈと乾度χは、次の式、ｈ＝ｆ（ｐ，χ）及びχ＝ｆ（ｐ，ｈ）から逆算する必要がある。以下の式（２４）、式（２５）は、それぞれ非共沸冷媒Ｒ４１０Ａ、Ｒ４０７Ｃのための演算式である。

ここで、Ｒ４１０Ａの場合ｐ_r＝ｐ_c／ｐとなり、Ｒ４０７Ｃの場合ｐ_r＝１０００ｐ_c／ｐとなる。但し、ｐ_cは臨界圧力、ｐは所与の二相圧力である。

以下の式において、エンタルピｈ、エントロピｓ及び温度Ｔは、次の式、ｈ＝ｆ（ｐ，χ）、ｓ＝ｆ（ｐ，ｈ）及びＴ＝ｆ（ｐ，ｈ）から逆算する必要がある。以下の式は、それぞれ非共沸冷媒Ｒ４１０Ａ、Ｒ４０７Ｃで共用の演算式である。

以上の式（２４）、式（２５）、式（２７）及び式（２８）によって、エンタルピｈ、乾度χ、エントロピｓ及び温度Ｔを得ることができる。

［サブクール領域の冷媒温度特性を求める場合］
サブクール冷媒温度特性の特徴は、飽和液冷媒温度特性のそれと同じであるので、サブクール冷媒温度特性を求めるための演算式は、間接的に飽和液冷媒温度特性を計算するための方程式を用いることで構成できる。なお、ここで示す演算方法は、冷媒が純冷媒であっても混合冷媒であっても使用することができる演算方法である。
一般的な実施例では、サブクール状態における圧力と温度は既知で、エンタルピｈ_scの計算が必要である。従って、サブクール状態におけるエンタルピｈ_scを既知のサブクール圧力ｐ_scとサブクール温度Ｔ_scによって計算する。そのための計算式として、次式を用いる。
ｈ_sc＝ｈ_ls−Ｃp_ls（Ｔ_ls−Ｔ_sc） （30）

この式（３０）を演算する流れを図１２のフローチャートに示す。この図１２において、先ず、（S1201）において既知のサブクール圧力ｐ_scとサブクール温度Ｔ_scを入力する。次に、（S1202）において既知のサブクール圧力ｐ_scを用いて飽和液温度Ｔ_lsを求める。飽和液温度Ｔ_lsが求まったら、これを用いて（S1203）において、飽和液エンタルピｈ_lsと飽和液定圧比熱Ｃp_lsを求める。これら値が求まったら、（S1204）において、式（３０）を用いて演算を行い、（S1205）でサブクールエンタルピｈ_scを得る。なお、図１２中のＴ_ls＝ｆ（ｐ_sc）、ｈ_ls＝ｆ（Ｔ_ls）及びＣp_ls＝ｆ（Ｔ_ls）は、それぞれのパラメータを求めるための高速演算式であり、前記実施例に示した方法で導出した陽方程式を用いるものである。

このように、それぞれ冷媒の種類や冷媒の状態ごとに異なる回帰式を適用することで、最適な高速演算式としての陽方程式を得ることが出来た。これらの方法で得た高速演算式を用いた計算は、従来の陰方程式を反復法で解く計算方法に比べて非常に計算速度の速いものとなっている。従来の計算方法の代表的なものとして、NISTのREFPROP6.01が存在するが、このREFPROPを用いた計算と本発明の高速演算式を用いた計算方法との計算時間の比較を行った。図１４〜図１６は、計算時間の比較結果を示したものである。なお、図１４において使用しているサブクール領域及び飽和領域における冷媒熱物性値を計算するための関数の定義を図１７に、同様に、図１５において使用しているスーパーヒート領域における冷媒熱物性値を計算するための関数の定義を図１８に示す。

図１４は、サブクール領域及び飽和領域における様々な冷媒熱物性値について計算時間の比較を行ったもので、REFPROPによる計算時間をｔ₁（ｓ）とし、本発明の高速演算式による計算時間をｔ₂（ｓ）としたときに、ｔ₁／ｔ₂によって計算時間の比率（倍率）を求めている。この図１４からも分かるように、全体としての計算時間は、従来の計算方法に比べて本発明の高速演算式による計算時間の方が、１７１．３倍速いことが分かる。同様に、図１５は、スーパーヒート領域における様々な冷媒熱物性値について計算時間の比較を行ったもので、この場合には、従来の計算方法に比べて本発明の高速演算式による計算時間の方が、４５．２倍速いことが分かる。さらに、図１６は、混合冷媒の場合の様々な冷媒熱物性値について計算時間の比較を行ったもので、この場合には、従来の計算方法に比べて本発明の高速演算式による計算時間の方が、Ｒ４１０Ａについては１４４９倍速く、Ｒ４０７Ｃについては１７９３倍速いことが分かる。
このように、従来の計算方法に比べて本発明の高速演算式による計算時間の方が、圧倒的に計算速度が速く、冷媒熱物性値の計算を行うシミュレータとしての実用性は圧倒的に高いことが分かる。

また、上記実施例では、フローチャートに基づいて説明を行ってきたが、本発明はこれに限られるものではなく、冷媒熱物性値の計算を行う装置としても構成することが出来る。例えば、図１３に示すように、全体を、ソースデータ演算処理部１１、回帰処理部１２、解式及び根選択処理部１３、冷媒熱物性演算部１４、及び、メモリ部１５で構成することが出来る。先ず、冷媒の指定、求めたい冷媒熱物性値の指定、及び、温度範囲とデータ点数の指定が、ソースデータ演算処理部１１に対して行われ、このソースデータ演算処理部１１でソースデータを演算によって求め、この求めたソースデータをメモリ部１５のソースデータ記憶部１６に記憶する。次に、回帰処理部１２において、求めたい冷媒熱物性値のソースデータをソースデータ記憶部１６から読出し、回帰式に当てはめて回帰処理を行い、回帰式の係数を得る。ここで得た複数の係数は、一時的に、演算過程一時記憶部１７に記憶される。そして、解式及び根選択処理部１３では、演算過程一時記憶部１７に記憶された係数を読み出して回帰式に当てはめ、この回帰式を解いて複数の根を得て、これらの根の中から有効な根を選択して、陽方程式記憶部１８に記憶する。冷媒熱物性演算部１４では、この陽方程式記憶部１８に記憶された陽方程式を高速演算式として用いて求めたい冷媒熱物性値の計算を行う。このように、冷媒熱物性値の計算を行う装置を構成することも出来る。また、このような装置においては、例えば、ソースデータ記憶部１６をデータベースとすることで、予めソースデータを有するような構成とすることもできる。

前記実施例１〜５により、冷媒の種類や状態に合わせてそれぞれ異なる演算方法で陽方程式を導出することが可能となったが、前述の通り陽方程式を導出するためには、求めたいパラメータを適宜変形した後に、式（３）等の３次方程式を用いて回帰処理を行って係数を得て、その係数を３次方程式に代入して解いて、その根の中から有効な根を選択して陽方程式を得ていた。このように、実際は求めたいパラメータ毎にそれぞれ異なる処理を行って、異なる回帰式によって回帰処理を行うものであり、この実施例においては、純冷媒Ｒ２２、Ｒ１３４ａと、混合冷媒Ｒ４１０Ａ、Ｒ４０７Ｃとについて、それぞれのパラメータ毎の陽方程式を導出するための詳細な演算式を示すとともに、これらを一覧表示して対応関係の把握を容易にした。なお、この実施例における演算式に適用する温度範囲については、Ｒ２２、Ｒ１３４ａ及びＲ４０７Ｃに関する飽和温度は−６０℃〜＋８０℃であり、Ｒ４１０Ａに関する飽和温度範囲は−６０℃〜＋６０℃であり、また、Ｒ２２、Ｒ１３４ａ、Ｒ４１０Ａ及びＲ４０７Ｃに関するスーパーヒート領域は０℃〜＋６５℃である。ただし、Ｒ４１０ＡのＴ_vs＝ｆ（ｈ）に関する−６０℃から＋２９℃の飽和温度と、Ｒ４０７ＣのＴ_vs＝ｆ（ｈ）に関する−５３℃から＋８０℃の飽和温度については除くものとする。

［純冷媒（Ｒ２２、Ｒ１３４ａ）が飽和状態の場合］
Ｒ２２、Ｒ１３４ａ等の純冷媒が飽和状態の場合における各パラメータの関係を表す陰方程式は、前記実施例１で述べたように、式（３）によって表現することが可能となる。よって、この式（３）を用いて回帰処理を行って陽方程式を導出するが、その際に、この式（３）におけるｕとｖに代入する各パラメータを適切な形に変形する必要がある。図１９に示すのは、Ｒ２２とＲ１３４ａが飽和状態の場合の主要パラメータの変形例であり、Ｒ２２に関するパラメータ毎の陰方程式No.をＩ１〜Ｉ９とし、Ｒ１３４ａに関するパラメータ毎の陰方程式No.をＩ１０〜Ｉ１８とした時に、各陰方程式に適用するパラメータの変形例を示している。例えば、Ｒ２２の飽和圧力ｐ_sと飽和温度Ｔ_sとの関係を表した陰方程式をＩ１とすると、式（３）の３次方程式のｕとｖにそれぞれ(ｐ_s−37.5049)／100とＴ_s−213.150とを代入したものが陰方程式Ｉ１となり、この陰方程式Ｉ１に回帰処理を行って各係数を得る。この図１９からも分かるように、パラメータ毎にそれぞれ異なる変形処理を行った後に回帰処理を行っている。

上記のように、最適な形式のパラメータを代入した陰方程式に対して回帰処理を行うことで、その陰方程式の係数を得ることができる。図２０は、Ｒ２２が飽和状態の場合の各陰方程式Ｉ１〜Ｉ９に回帰処理を行って求めた各係数を表したものであり、また、図２１は、Ｒ１３４ａが飽和状態の場合の各陰方程式Ｉ１０〜Ｉ１８に回帰処理を行って求めた各係数を表したものである。このようにして求めた係数を用いて陽方程式を導出する。図２２は、Ｒ２２が飽和状態の場合の各陰方程式Ｉ１〜Ｉ９から導出した陽方程式と、Ｒ１３４ａが飽和状態の場合の各陰方程式Ｉ１０〜Ｉ１８から導出した陽方程式とを示したものであり、また、これらの陽方程式を導出するのに用いた各式を示している。

例えば、Ｒ２２が飽和状態の場合の飽和圧力ｐ_sと飽和温度Ｔ_sとの関係を表した陽方程式を導出するには、陰方程式Ｉ１を回帰処理して係数ａ_１〜ａ_８を求める。ここで、陰方程式（３）は式（４）のように整理でき、この式（４）における各係数を図９（Ａ−６）に示すａ〜ｄで置換え（但し、式中のｘ＝ｖ）、ｕ＝ｙで置換えると、式（４）は以下のようになる。
ａｙ^３＋ｂｙ^２＋ｃｙ＋ｄ＝０ （３１）
この式（３１）の根のうち有効な根となるのが図８（Ａ−２）であり、これによって有効な根ｙを求める。このｙを図２２に示すＲ２２の場合のｐ_sとＴ_sとの関係を表した陽方程式に代入することで、最終的な陽方程式を得る。このように、図２２に示した数式を用いることで、Ｒ２２が飽和状態の場合の各パラメータに関する陰方程式Ｉ１〜Ｉ９からそれぞれ陽方程式を導出することができ、同様に、Ｒ１３４ａが飽和状態の場合の各パラメータに関する陰方程式Ｉ１〜Ｉ９からそれぞれ陽方程式を導出することができる。

［混合冷媒（Ｒ４１０Ａ、Ｒ４０７Ｃ）が飽和状態の場合］
混合冷媒であるＲ４１０ＡとＲ４０７Ｃが飽和状態の場合についても同様に、図２３に示すのは、Ｒ４１０ＡとＲ４０７Ｃが飽和状態の場合における主要パラメータの変形例であり、Ｒ４１０Ａに関するパラメータ毎の陰方程式No.をＩ１９〜Ｉ２８とし、Ｒ４０７Ｃに関するパラメータ毎の陰方程式No.をＩ２９〜Ｉ３８とした時に、各陰方程式に適用するパラメータの変形例を示している。Ｒ４１０Ａに関する陰方程式Ｉ１９〜Ｉ２８をそれぞれ回帰処理して求めた係数が図２４であり、Ｒ４０７Ｃに関する陰方程式Ｉ２９〜Ｉ３８をそれぞれ回帰処理して求めた係数が図２５である。これらに基づいて導出したＲ４１０ＡとＲ４０７Ｃが飽和状態の場合の陽方程式を図２６に示し、また、これらの陽方程式を導出するのに用いた各式を図２６に示す。

［冷媒がスーパーヒート状態の場合］
冷媒がスーパーヒート状態の場合における各パラメータの関係を表す陰方程式は、前記実施例２で述べたように、式（１８）によって表現することが可能となる。この式（１８）のパラメータの具体的な変形パターンとして、図２７（ａ）は、Ｒ２２とＲ１３４ａがスーパーヒート状態の場合における主要パラメータの変形例を示したものであり、Ｒ２２に関するパラメータ毎の陰方程式No.をＩ３９〜Ｉ４１とし、Ｒ１３４ａに関するパラメータ毎の陰方程式No.をＩ４２〜Ｉ４４としている。また、図２７（ｂ）は、Ｒ４１０ＡとＲ４０７Ｃがスーパーヒート状態の場合における主要パラメータの変形例を示したものであり、Ｒ４１０Ａに関するパラメータ毎の陰方程式No.をＩ４５〜Ｉ４７とし、Ｒ４０７Ｃに関するパラメータ毎の陰方程式No.をＩ４８〜Ｉ５０としている。

上記Ｒ２２及びＲ１３４ａに関する陰方程式Ｉ３９〜Ｉ４４をそれぞれ回帰処理して求めた係数が図２８であり、Ｒ４１０Ａ及びＲ４０７Ｃに関する陰方程式Ｉ４５〜Ｉ５０をそれぞれ回帰処理して求めた係数が図２９である。これらに基づいて導出したＲ２２及びＲ１３４ａがスーパーヒート状態の場合の陽方程式及び陽方程式を導出するのに用いた各式を図３０（ａ）に示し、Ｒ４１０Ａ及びＲ４０７Ｃがスーパーヒート状態の場合の陽方程式及び陽方程式を導出するのに用いた各式を図３０（ｂ）に示す。

［混合冷媒（Ｒ４１０Ａ、Ｒ４０７Ｃ）の二相温度特性を求める場合］
混合冷媒の二相温度特性を求める場合における各パラメータの関係を表す陰方程式は、図３３に示す式（３２）によって表現することが可能となる。この式（３２）におけるＡは、前記式（２９）によって与えられるものであり、この式（３２）のパラメータであるｕとｖ、及び、式（２９）におけるパラメータであるｐ_rの具体的な変形パターンを示したものが図３１である。この図３１において、Ｒ４１０Ａに関するパラメータ毎の陰方程式No.をＩ５１〜Ｉ５３とし、Ｒ４０７Ｃに関するパラメータ毎の陰方程式No.をＩ５４〜Ｉ５６としている。

上記Ｒ４１０Ａ及びＲ４０７Ｃに関する陰方程式Ｉ５１〜Ｉ５６をそれぞれ回帰処理して求めた係数が図３２である。これらに基づいて導出した混合冷媒の二相温度特性を求める場合の陽方程式を図３３に示す。なお、陽方程式の導出のために、図３３の式（３２）又は式（３３）を用いている。

このように、純冷媒Ｒ２２、Ｒ１３４ａと、混合冷媒Ｒ４１０Ａ、Ｒ４０７Ｃとについて、それぞれの陰関数毎にパラメータを最適な形式に変形して回帰処理を行い、それによって得た係数を用いて具体的な陽方程式を導出する過程について説明したが、この実施例６において導出したＲ２２、Ｒ１３４ａ、Ｒ４１０Ａ及びＲ４０７Ｃに関する陽方程式の計算精度及び計算速度について図面を用いて説明する。

図３４に示すのは、Ｒ２２、Ｒ１３４ａ、Ｒ４１０Ａ及びＲ４０７Ｃに関する陽方程式の計算精度を検証するために、NISTのREFPROP6.01によって計算した値と、この実施例６に示した陽方程式によって計算した値とを比較し、その場合の平均誤差と最大誤差とを各陽方程式ごとに求めた。飽和状態の陽方程式については140点を計算して求め、スーパーヒート状態の陽方程式については3000点を計算して求め、二相温度特性の陽方程式については4000点を計算して求めた。この図３４からも分かるように、各冷媒の高速演算式の総平均相対誤差は0.021％以下であることがわかる。またρ_sp =ｆ（ｐ，Ｔ）の最大相対誤差は2.0％以下と大きくなっているが、他の高速演算式の最大相対誤差は0.61％以下である。

また、図３５に示すのは、Ｒ２２、Ｒ１３４ａ、Ｒ４１０Ａ及びＲ４０７Ｃに関する陽方程式の計算速度を検証するために、陽方程式ごとに、NISTのREFPROP6.01による計算とこの実施例６に示した陽方程式による計算とのそれぞれについて平均計算時間と最大計算時間とを求めるとともに、REFPROPによる計算に比べてこの実施例６に示した陽方程式による計算がどの程度高速化されたかを示す倍率を求めた。その結果、図３５からも分かるように、純冷媒Ｒ２２及びＲ１３４ａに関する平均計算速度は、REFPROPの速度と比較してそれぞれ約140倍高速になっており、混合冷媒Ｒ４１０Ａ及びＲ４０７Ｃに関する平均計算速度は、REFPROPの速度と比較してそれぞれ約1000倍高速になっていることが分かる。

本発明の冷媒熱物性値の計算において使用する陽方程式を導出する流れを示したフローチャート図である。 図１の（S103）における処理の詳細を示したフローチャート図である。 純冷媒の飽和領域における冷媒熱物性を求める場合の図１の（S105）における処理の詳細を示したフローチャート図である。 純冷媒の飽和領域における冷媒熱物性を求める場合の図１の（S111）における処理の詳細を示したフローチャート図である。 純冷媒の飽和領域における冷媒熱物性を求める場合の図４の（S402）における処理の詳細を示したフローチャート図である。 純冷媒の飽和領域における冷媒熱物性を求める場合の図４の（S403）における処理の詳細を示したフローチャート図である。 純冷媒の飽和領域における冷媒熱物性を求める場合の図４の（S405）における処理において図５を引用する場合に、その図５における（S503）の処理に置換えて使用するステップ図である。 図５の（S506）における処理で使用する演算式である。 図５の（S505）及び（S506）における処理で使用する演算式である。 純冷媒のスーパーヒート領域における冷媒熱物性を求める場合の図１の処理に置換えて使用するステップ図である。 純冷媒のスーパーヒート領域における冷媒熱物性を求める場合の図１の（S105）における処理の詳細を示したフローチャート図である。 サブクール領域の冷媒温度特性を求める場合の演算の流れを示したフローチャート図である。 図１における計算方法を実現するための計算装置の構成を示したブロック図である。 サブクール領域及び飽和領域における様々な冷媒熱物性値について、REFPROP6.01による計算時間と本発明の高速演算式による計算時間の比較を行った表である。 スーパーヒート領域における様々な冷媒熱物性値について、REFPROP6.01による計算時間と本発明の高速演算式による計算時間の比較を行った表である。 混合冷媒の場合の様々な冷媒熱物性値について、REFPROP6.01による計算時間と本発明の高速演算式による計算時間の比較を行った表である。 図１４において示したサブクール領域及び飽和領域における冷媒熱物性値を計算するための関数の定義を示した表である。 図１５において示したスーパーヒート領域における冷媒熱物性値を計算するための関数の定義を示した表である。 飽和冷媒Ｒ２２及びＲ１３４ａに関する陰方程式の主要パラメータの変形例を示した表である。 飽和冷媒Ｒ２２に関する陰関数方程式の係数を示した表である。 飽和冷媒Ｒ１３４ａに関する陰関数方程式の係数を示した表である。 飽和冷媒Ｒ２２及びＲ１３４ａに関する陽方程式を示した表である。 飽和冷媒Ｒ４１０Ａ及びＲ４０７Ｃに関する陰方程式の主要パラメータの変形例を示した表である。 飽和冷媒Ｒ４１０Ａに関する陰方程式の係数を示した表である。 飽和冷媒Ｒ４０７Ｃに関する陰方程式の係数を示した表である。 飽和冷媒Ｒ４１０Ａ及びＲ４０７Ｃに関する陽方程式を示した表である。 スーパーヒート冷媒に関する陰関数方程式の主要パラメータの変形例を示した表であり、（ａ）は純冷媒の場合、（ｂ）は混合冷媒の場合を表している。 純冷媒がスーパーヒート状態の場合の陰方程式の係数を示した表である。 混合冷媒がスーパーヒート状態の場合の陰方程式の係数を示した表である。 スーパーヒート冷媒に関する陽方程式を示した表であり、（ａ）は純冷媒の場合、（ｂ）は混合冷媒の場合を表している。 二相冷媒Ｒ４１０Ａ及びＲ４０７Ｃに関する陰方程式の主要パラメータの変形例を示した表である。 二相冷媒Ｒ４１０Ａ及びＲ４０７Ｃに関する陰方程式の係数を示した表である。 二相冷媒Ｒ４１０Ａ及びＲ４０７Ｃに関する陽方程式を示した表、及び、この陽方程式を導出するのに用いた数式を表したものである。 Ｒ２２、Ｒ１３４ａ、Ｒ４１０Ａ及びＲ４０７Ｃについて、NISTのREFPROP6.01によって計算した値と、実施例６に示した陽方程式によって計算した値とを比較した場合の相対誤差を示した表である。 Ｒ２２、Ｒ１３４ａ、Ｒ４１０Ａ及びＲ４０７Ｃについて、NISTのREFPROP6.01による計算と、実施例６に示した陽方程式による計算とのそれぞれの計算速度を示すとともに、両者の計算速度の比較を行った結果を示した表である。

符号の説明

（S101）〜（S116）…図１における各ステップ、（S201）〜（S206）…図２における各ステップ、（S301）〜（S310）…図３における各ステップ、（S401）〜（S407）…図４における各ステップ、（S501）〜（S511）…図５における各ステップ、（S601）〜（S608）…図６における各ステップ、（S701）…図７におけるステップ、（S1001）〜（S1002）…図１０における各ステップ、（S1101）〜（S1111）…図１１における各ステップ、（S1201）〜（S1205）…図１２における各ステップ、１１…ソースデータ演算処理部、１２…回帰処理部、１３…解式及び根選択処理部、１４…冷媒熱物性演算部、１５…メモリ部、１６…ソースデータ記憶部、１７…演算過程一時記憶部、１８…陽方程式記憶部。

Claims (12)

1. 演算処理部及びメモリ部を有したコンピュータにおいて、冷媒熱物性値のソースデータを演算処理部で演算して得てメモリ部に記憶するソースデータ演算処理手段と、特定の高次多項式をメモリ部に記憶された前記冷媒熱物性値のソースデータのうち、求めたい冷媒熱物性値によって演算処理部における演算により回帰処理することで前記高次多項式の係数を得てメモリ部に記憶する回帰処理手段と、前記メモリ部に記憶された係数を代入した高次多項式を演算処理部における演算により解くことで前記高次多項式の根を得てメモリ部に記憶し、さらに前記根の中から有効な根を選択してメモリ部に記憶する解式及び根選択処理手段とによって、求めたい冷媒熱物性値の関係を表した陽方程式を得て、この陽方程式を用いて冷媒熱物性値の計算を演算処理部において行うことを特徴とする冷媒熱物性値の計算方法。
2. ソースデータ演算処理手段において求めるソースデータは、冷媒熱物性値の関係を表した陰方程式を演算処理部における演算により反復法によって解いて求めたものを含むことを特徴とする請求項１記載の冷媒熱物性値の計算方法。
3. ソースデータ演算処理手段では、指定した冷媒の温度範囲について演算処理部における演算によりソースデータを求め、回帰処理手段では、前記冷媒の温度範囲を指定して求めたソースデータのうち、求めたい冷媒熱物性値について高次多項式に当てはめて最小二乗法により演算処理部における演算により回帰処理し、この回帰処理が収束した場合に前記高次多項式の係数を得てメモリ部に記憶するようにしたことを特徴とする請求項１又は２記載の冷媒熱物性値の計算方法。
4. 回帰処理手段において行った回帰処理が収束しなかった場合には、冷媒の温度範囲を調整し、その調整後の温度範囲のソースデータを用いて演算処理部における演算により回帰処理を行うようにし、温度範囲を調整しながら回帰処理を設定した回数行っても収束しない場合には、回帰処理に用いたソースデータに対して変換式を適用した後に再度演算処理部における演算により回帰処理を行うようにしたことを特徴とする請求項３記載の冷媒熱物性値の計算方法。
5. 解式及び根選択処理手段では、求めた前記高次多項式の根をそれぞれソースデータと演算処理部において比較して、最も誤差の小さい根を有効な根としてメモリ部に記憶することを特徴とする請求項１記載の冷媒熱物性値の計算方法。
6. 請求項１乃至５記載のソースデータ演算処理手段、回帰処理手段、及び、解式及び根選択処理手段とをコンピュータに実行させ、これによって求めたい冷媒熱物性値の関係を表した陽方程式を得て、この陽方程式を用いて冷媒熱物性値の計算を行うことを特徴とする冷媒熱物性値の計算プログラム。
7. 請求項１乃至５記載のソースデータ演算処理手段、回帰処理手段、及び、解式及び根選択処理手段とをコンピュータに実行させ、これによって求めたい冷媒熱物性値の関係を表した陽方程式を得て、この陽方程式を用いて冷媒熱物性値の計算を行うことを特徴とする冷媒熱物性値の計算プログラムを記憶した記憶媒体。
8. 冷媒熱物性値のソースデータを演算処理によって得てメモリ部に記憶するソースデータ演算処理部と、特定の高次多項式を前記メモリ部に記憶された前記冷媒熱物性値のソースデータのうち、求めたい冷媒熱物性値によって演算により回帰処理することで前記高次多項式の係数を得てメモリ部に記憶する回帰処理部と、前記メモリ部に記憶された係数を代入した高次多項式を解くことで前記高次多項式の根を得てメモリ部に記憶し、さらに前記根の中から有効な根を選択してメモリ部に記憶する解式及び根選択処理部とを具備し、これらの処理部を経ることで求めたい冷媒熱物性値の関係を表した陽方程式を得てこれをメモリ部に記憶し、このメモリ部に記憶された陽方程式を用いて冷媒熱物性値の計算を行う冷媒熱物性演算部を具備したことを特徴とする冷媒熱物性値の計算装置。
9. ソースデータ演算処理部において求めるソースデータは、冷媒熱物性値の関係を表した陰方程式を演算により反復法によって解いて求めたものを含むことを特徴とする請求項８記載の冷媒熱物性値の計算装置。
10. ソースデータ演算処理部では、指定した冷媒の温度範囲について演算によりソースデータを求め、回帰処理部では、前記冷媒の温度範囲を指定して求めたソースデータのうち、求めたい冷媒熱物性値について高次多項式に当てはめて最小二乗法により演算により回帰処理し、この回帰処理が収束した場合に前記高次多項式の係数を得てメモリ部に記憶するようにしたことを特徴とする請求項８又は９記載の冷媒熱物性値の計算装置。
11. 回帰処理部において行った回帰処理が収束しなかった場合には、冷媒の温度範囲を調整し、その調整後の温度範囲のソースデータを用いて演算により回帰処理を行うようにし、温度範囲を調整しながら回帰処理を設定した回数行っても収束しない場合には、回帰処理に用いたソースデータに対して変換式を適用した後に再度演算により回帰処理を行うようにしたことを特徴とする請求項１０記載の冷媒熱物性値の計算装置。
12. 解式及び根選択処理部では、求めた前記高次多項式の根をそれぞれソースデータと比較して、最も誤差の小さい根を有効な根としてメモリ部に記憶することを特徴とする請求項８記載の冷媒熱物性値の計算装置。

Images