(57)【要約】      （修正有）【課題】プラントの流路網の計算について、プラント構成が一部変更された場合においても変更の手間を最小限とする。【解決手段】シミュレーション装置１００は、入力部３１と、流路系構成記憶部２１と、流量バランスベクトル設定部１１と、流量バランスベクトル記憶部２２と、変分比演算部１２と、変分比記憶部２３と、収束演算により複数の圧力および複数の流量の値を算出するとともに収束判定を行う収束演算部１３と、結果を記憶する収束状態記憶部２４と、演算結果を出力する出力部３２とを有する。流量バランスベクトル設定部１１は、流量バランス式にもとづき流量バランスベクトルを設定する。変分比演算部１２は、複数の圧力の所定の微小変動に対する流量の変動分のベクトルを算出する。【選択図】図１
2017041124.tif 000002
【特許請求の範囲】
【請求項１】
解析対象とするプラントの流路系についての複数の評価点に対応する複数の圧力および複数の流量の値を算出するシミュレーション装置であって、
前記流路系の構成、各部のパラメータおよび境界条件に関するデータを入力として外部から受け入れる入力部と、
前記入力部が受け入れた前記流路系の構成に関するデータに基づいて、前記流路系の構成を記憶する流路系構成記憶部と、
前記流路系構成記憶部に記憶された前記流路系の構成に基づいて前記複数の圧力の関数として表される前記複数の流量について、当該複数の流量間のバランス式にもとづいて流量バランスベクトルを設定する流量バランスベクトル設定部と、
前記流量バランスベクトル設定部で設定された流量バランスベクトルを記憶する流量バランスベクトル記憶部と、
前記流量バランスベクトル記憶部に記憶された前記流量バランスベクトルを用いて、変分比を演算する変分比演算部と、
前記変分比演算部で算出された変分比を記憶する変分比記憶部と、
前記変分比記憶部に記憶された変分比を用いて、前記複数の圧力の圧力推定値を設定して評価値を算出し、前記評価値に基づいて前記圧力推定値を順次更新する収束演算により、前記複数の圧力の値を算出し、当該複数の圧力から前記複数の流量の値を算出するとともに収束判定を行う収束演算部と、
前記収束演算部の各サイクルでの演算結果を更新しながら記憶する収束状態記憶部と、
前記収束演算部が収束と判定した場合に、前記収束状態記憶部に記憶する演算結果を出力する出力部と、
を備え、
前記流量バランスベクトル設定部における設定は次のバランス式（１）を設定することにより行い、
Ｆ（ｐ）＝０            （１）
前記変分比演算部における前記変分比の演算は、次の式（２）に表す、前記複数の圧力の所定の微小変動に対する前記流量の変動分のベクトルを算出することにより行うことを特徴とするシミュレーション装置。
Ａ（ｘ）ｐ＝［Ｆ（ｘ＋σｐ）−Ｆ（ｘ）］／σ            …（２）
ただし、式（１）、（２）で、ｐは、前記複数の圧力を要素とする圧力ベクトル、Ｆ（ｐ）は、前記圧力ベクトルの各要素の関数として表される複数の流量バランス式を要素とする流量バランスベクトル、０は、ゼロベクトルを表し、
また、式（２）で、Ａ（ｘ）は、ベクトルＦ（ｐ）の各要素のベクトルｐの各要素による偏微分を要素とするヤコビアン行列、σは所定の微小量を表す。
【請求項２】
前記流路系構成記憶部は、前記流量バランスベクトルの設定のための、前記流路系のそれぞれのノード点について流入項目および流出項目を対応させた流量バランステーブルを有することを特徴とする請求項１に記載のシミュレーション装置。
【請求項３】
前記流量バランスベクトル設定部は、前記流量バランスベクトルの各流量バランス要素の数値のオーダーのばらつきを抑制するために、それぞれの要素が有する流量係数に基づいて設定される係数を各流量バランス要素に乗ずる乗算部を有することを特徴とする請求項１または請求項２に記載のシミュレーション装置。
【請求項４】
解析対象とするプラントの流路系についての複数の評価点での複数の圧力および複数の流量を算出するシミュレーション方法であって、
入力部が、前記流路系の構成、各部のパラメータおよび境界条件に関するデータを入力として外部から受け入れる入力ステップと、
流路系構成記憶部が、前記入力部が受け入れた前記流路系の構成に関するデータに基づいて、前記流路系の構成を記憶する流路系構成記憶ステップと、
流量バランスベクトル設定部が、前記流路系構成記憶部に記憶された前記流路系の構成に基づいて前記複数の圧力の関数として表される前記複数の流量について、当該複数の流量間のバランス式にもとづいて流量バランスベクトルを設定する流量バランスベクトル設定ステップと、
流量バランスベクトル記憶部が、前記流量バランスベクトル設定ステップで設定された流量バランスベクトルを記憶する流量バランスベクトル記憶ステップと、
変分比演算部が、前記流量バランスベクトル記憶部に記憶された前記流量バランスベクトルを用いて、変分比を演算する変分比演算ステップと、
変分比記憶部が、前記変分比演算ステップで算出された変分比を記憶する変分比記憶ステップと、
収束演算部が、前記変分比記憶部に記憶された変分比を用いて、前記複数の圧力の圧力推定値を設定して評価値を算出し、前記評価値に基づいて前記圧力推定値を順次更新する収束演算により、前記複数の圧力の値を算出し、当該複数の圧力から前記複数の流量の値を算出するとともに収束判定を行う収束演算ステップと、
前記収束演算ステップに並行して、収束状態記憶部が、前記収束演算部の各サイクルでの演算結果を更新しながら記憶する収束状態記憶ステップと、
前記収束演算部が収束と判定した場合に、出力部が、前記収束状態記憶部に記憶する演算結果を出力する出力ステップと、
を有することを特徴とするシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、プラント等の流路系のシミュレーション装置およびシミュレーション方法に関する。
【背景技術】
【０００２】
発電プラントや工場設備など（以下、プラントという。）の運転操作の訓練、および発生した事象の解析のためには、プラント内の流路網を流れるガスや流体などの流量、圧力およびエンタルピーなどを求める必要がある。その際、そのシミュレーションモデルは、プラント内の系統や工程ごとに作成される。プラントの構成上、系統間の連結点は少なく、また、事象を限定した場合は、独立に作られる場合が多い。
【０００３】
一方、事故事象に関しては、たとえば、新たに破断箇所を考慮する必要性が出るなどプラントモデルの修正が必要となる場合が多い。流路網の計算では、一般的には圧力と流量に関する非線形の方程式の集まりをニュートン法で解く場合が多い。この際、圧力点での流量に関する圧力の微係数行列（ヤコビアン）の逆行列を求めるなど、大きなサイズの行列を扱う必要がある。また、連結点が少ない場合は、この行列は疎なものになり、ゼロの要素が多いが、大きなメモリーサイズを必要とすることになる。その場合、モデルの変更に伴う処理系の変更は複雑なものになり、行列が疎である割に多大な手間を要するようになる。
【０００４】
流路網の計算としては、特許文献１には換気系の計算においてボリュームノード間の流量を、等価流路網を設定して計算する方法が示されている。換気系統において高速化する方法が示されている。
【０００５】
非特許文献１には管路系のノードジャンクション法でのシミュレーション例が示されている。流路網（管路網）の計算には反復法あるいはニュートン法が使われていることが示されている。
【０００６】
ニュートン法で逆行列を陽に計算しない方法として、非特許文献２にはＪＦＮＫ法（Ｊａｃｏｉａｎ  Ｆｒｅｅ  Ｎｅｗｔｏｎ−Ｋｒｙｌｏｖ  Ｍｅｔｈｏｄ）が提案、応用されている。ＪＦＮＫ法は、クリロフ法に基づく逆行列計算の際にはヤコビアンを直接計算しなくて済むことから、メモリーを多く使わない方法として提案されたものである。しかしながら、これらはモデルの変更に伴う処理系の変更を容易にすることを対象としたものではなく、また、流路網を対象にしたものでもない。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特許第４８９８３０６号公報
【非特許文献】
【０００８】
【非特許文献１】富士総合研究所編「管路内の流れのシミュレーションプログラム」１９９５年９月
【非特許文献２】日本機械学会論文（Ｂ編）７５巻７５９号「完全保存形差分スキームとＪＦＮＫ法による非圧縮性流れの非分離解法」
【発明の概要】
【発明が解決しようとする課題】
【０００９】
本発明の実施形態はこのような事情によりなされたもので、プラントの流路網の計算について、プラント構成が一部変更された場合においても、シミュレーションにおける変更の手間を最小限とすることを目的とする。
【課題を解決するための手段】
【００１０】
上述の目的を達成するため、本実施形態は、解析対象とするプラントの流路系についての複数の評価点に対応する複数の圧力および複数の流量の値を算出するシミュレーション装置であって、前記流路系の構成、各部のパラメータおよび境界条件に関するデータを入力として外部から受け入れる入力部と、前記入力部が受け入れた前記流路系の構成に関するデータに基づいて、前記流路系の構成を記憶する流路系構成記憶部と、前記流路系構成記憶部に記憶された前記流路系の構成に基づいて前記複数の圧力の関数として表される前記複数の流量について、当該複数の流量間のバランス式にもとづいて流量バランスベクトルを設定する流量バランスベクトル設定部と、前記流量バランスベクトル設定部で設定された流量バランスベクトルを記憶する流量バランスベクトル記憶部と、前記流量バランスベクトル記憶部に記憶された前記流量バランスベクトルを用いて、変分比を演算する変分比演算部と、前記変分比演算部で算出された変分比を記憶する変分比記憶部と、前記変分比記憶部に記憶された変分比を用いて、前記複数の圧力の圧力推定値を設定して評価値を算出し、前記評価値に基づいて前記圧力推定値を順次更新する収束演算により、前記複数の圧力の値を算出し、当該複数の圧力から前記複数の流量の値を算出するとともに収束判定を行う収束演算部と、前記収束演算部の各サイクルでの演算結果を更新しながら記憶する収束状態記憶部と、前記収束演算部が収束と判定した場合に、前記収束状態記憶部に記憶する演算結果を出力する出力部と、を備え、前記流量バランスベクトル設定部における設定はバランス式Ｆ（ｐ）＝０…（１）を設定することにより行い、前記変分比演算部における前記変分比の演算は、Ａ（ｘ）ｐ＝［Ｆ（ｘ＋σｐ）−Ｆ（ｘ）］／σ…（２）に表す、前記複数の圧力の所定の微小変動に対する前記流量の変動分のベクトルを算出することにより行うことを特徴とする。ただし、式（１）、（２）で、ｐは、前記複数の圧力を要素とする圧力ベクトル、Ｆ（ｐ）は、前記圧力ベクトルの各要素の関数として表される複数の流量バランス式を要素とする流量バランスベクトル、０は、ゼロベクトルを表し、また、式（２）で、Ａ（ｘ）は、ベクトルＦ（ｐ）の各要素のベクトルｐの各要素による偏微分を要素とするヤコビアン行列、σは所定の微小量を表す。
【００１１】
また、本実施形態は、解析対象とするプラントの流路系についての各評価点での複数の圧力および複数の流量を算出するシミュレーション方法であって、入力部が、前記流路系の構成、各部のパラメータおよび境界条件に関するデータを入力として外部から受け入れる入力ステップと、流路系構成記憶部が、前記入力部が受け入れた前記流路系の構成に関するデータに基づいて、前記流路系の構成を記憶する流路系構成記憶ステップと、流量バランスベクトル設定部が、前記流路系構成記憶部に記憶された前記流路系の構成に基づいて前記複数の圧力の関数として表される前記複数の流量について、当該複数の流量間のバランス式にもとづいて流量バランスベクトルを設定する流量バランスベクトル設定ステップと、流量バランスベクトル記憶部が、前記流量バランスベクトル設定ステップで設定された流量バランスベクトルを記憶する流量バランスベクトル記憶ステップと、変分比演算部が、前記流量バランスベクトル記憶部に記憶された前記流量バランスベクトルを用いて、変分比を演算する変分比演算ステップと、変分比記憶部が、前記変分比演算ステップで算出された変分比を記憶する変分比記憶ステップと、収束演算部が、前記変分比記憶部に記憶された変分比を用いて、前記複数の圧力の圧力推定値を設定して評価値を算出し、前記評価値に基づいて前記圧力推定値を順次更新する収束演算により、前記複数の圧力の値を算出し、当該複数の圧力から前記複数の流量の値を算出するとともに収束判定を行う収束演算ステップと、前記収束演算ステップに並行して、収束状態記憶部が、前記収束演算部の各サイクルでの演算結果を更新しながら記憶する収束状態記憶ステップと、前記収束演算部が収束と判定した場合に、出力部が、前記収束状態記憶部に記憶する演算結果を出力する出力ステップと、を有することを特徴とする。
【発明の効果】
【００１２】
本発明の実施形態によれば、プラントの流路網の計算について、プラント構成が一部変更された場合においても、シミュレーションにおける変更の手間を最小限とすることが可能となる。
【図面の簡単な説明】
【００１３】
【図１】実施形態に係るシミュレーション装置の構成を示すブロック図である。
【図２】プラントの流路系の構成例を示す接続関係図である。
【図３】実施形態に係るシミュレーション装置におけるプラントの流路系の構成例に関する流量バランステーブルである。
【図４】プラントの流路系の構成の例の結合後の構成を示す系統図である。
【図５】実施形態に係るシミュレーション装置における結合後の構成に関する流量バランステーブルである。
【図６】実施形態に係るシミュレーション装置の収束演算部における手順を示すフロー図である。
【図７】実施形態に係るシミュレーション方法の手順を示すフロー図である。
【図８】実施形態に係るシミュレーション方法の結合後の系統についての計算の手順を示すフロー図である。
【図９】従来のシミュレーション方法の手順を示すフロー図である。
【発明を実施するための形態】
【００１４】
以下、図面を参照して、本発明の実施形態に係るシミュレーション装置およびシミュレーション方法について説明する。
【００１５】
図１は、実施形態に係るシミュレーション装置の構成を示すブロック図である。シミュレーション装置１００は、解析対象とするプラントの流路系についての複数の評価点に対応する複数の圧力および複数の流量の値を算出する。シミュレーション装置１００は、演算部１０、記憶部２０、入力部３１および出力部３２を有する。
【００１６】
演算部１０は、流量バランスベクトル設定部１１、変分比演算部１２、収束演算部１３、および時間進展管理部１４を有する。また、記憶部２０は、流路系構成記憶部２１、流量バランスベクトル記憶部２２、変分比記憶部２３、および収束状態記憶部２４を有する。
【００１７】
入力部３１は、流路系の構成、各部のパラメータおよび境界条件に関するデータを入力として外部から受け入れる。流路系構成記憶部２１は、入力部３１が受け入れた流路系の構成、各部のパラメータおよび境界条件に関するデータを記憶する。
【００１８】
流量バランスベクトル設定部１１は、流路系構成記憶部２１に記憶された流路系の構成に基づいて複数の圧力の関数として表される複数の流量について、複数の流量間のバランス式にもとづいて流量バランスベクトルを設定する。
【００１９】
図２は、プラントの流路系の構成例を示す接続関係図である。図２で示されているのは、第１の流路と第２の流路の独立した２つの流路である。それぞれの流路は、両端が端点となっており、境界条件として与えられる。また、端点間には、バルブが配されており、それぞれの弁の前後が圧力ノードとなる。
【００２０】
すなわち、第１の流路は、端点Ｎ０で境界条件は圧力ｐ０、ノード点Ｎ１で圧力ｐ１、ノード点Ｎ２で圧力ｐ２、ノード点Ｎ３で圧力ｐ３、端点Ｎ４で境界条件は圧力ｐ４のように各ノード点での圧力が設定されている。このとき、端点Ｎ０とノード点Ｎ１間のバルブの流量はｗ０１、ノード点Ｎ１とノード点Ｎ２間のバルブの流量はｗ１２、ノード点Ｎ２とノード点Ｎ３間のバルブの流量はｗ２３、ノード点Ｎ３と端点Ｎ４間のバルブの流量はｗ３４である。
【００２１】
同様に、第２の流路は、端点Ｎ５で境界条件は圧力ｐ５、ノード点Ｎ６で圧力ｐ６、ノード点Ｎ７で圧力ｐ７、端点Ｎ８で境界条件は圧力ｐ８のように各ノード点での圧力が設定されている。このとき、端点Ｎ５とノード点Ｎ６間のバルブの流量はｗ５６、ノード点Ｎ６とノード点Ｎ７間のバルブの流量はｗ６７、ノード点Ｎ７と端点Ｎ８間のバルブの流量はｗ７８である。
【００２２】
図３は、実施形態に係るシミュレーション装置におけるプラントの流路系の構成例に関する流量バランステーブルである。流路系構成記憶部２１は、このテーブル形式で流路構成を記憶する。流量バランステーブルは、図３に示すように、各ノード点について、流入量、流出量をそれぞれ割り振っている。また、ノード点が端点の場合は、端点である旨が入力されている。また、変更の経緯を残すために経緯欄が設けられている。
【００２３】
たとえば、ノードＮ０は、端点であり、流入はなし、流出はｗ０１である。また、ノード点Ｎ１は、流入がｗ０１、流出がｗ１２である。
【００２４】
図４は、プラントの流路系の構成の例の結合後の構成を示す系統図である。すなわち、図２に示した構成との違いは、第２の流路の端点Ｎ５が端点でなくなくなり、第１の流路のノード点Ｎ１に結合された場合を示している。この結果、第１の流路と第２の流路は互いに独立ではなく、互いに結合された流路となっている。
【００２５】
図５は、結合後の構成に関する流量バランステーブルである。図５に示すように、端点Ｎ５の欄が削除される。また、端点Ｎ５から中間のノード点に変わったため、ノード点Ｎ１からノード点Ｎ６への流れｗ１６が生ずる。このため、ノード点Ｎ１の流出欄、ノード点Ｎ６の流入欄にそれぞれｗ１６が追加される。また、図２に示す構成から図４に示す構成に変わった経緯を、経緯欄に収納している。
【００２６】
流量バランスベクトル設定部１１は、それぞれの流路に関する流量バランスベクトルを設定する。流量バランスベクトルは、流量バランステーブルを用いて、各ノード点での流量のバランスに基づいて設定される。
【００２７】
たとえば、図２で示した構成における第１の流路に関するノード点Ｎ１、Ｎ２、およびＮ３のそれぞれにおいて流入流量から流出流量を減じるとゼロになる関係を次のように数式化する。
ノード点Ｎ１：    ｗ０１−ｗ１２＝０            …（１）
ノード点Ｎ２：    ｗ１２−ｗ２３＝０            …（２）
ノード点Ｎ３：    ｗ２３−ｗ３４＝０            …（３）
【００２８】
また、第２の流路に関するノード点Ｎ６、およびＮ７における流入流量から流出流量を減じるとゼロになる関係を次のように数式化する。
ノード点Ｎ６：    ｗ５６−ｗ６７＝０            …（４）
ノード点Ｎ７：    ｗ６７−ｗ７８＝０            …（５）
【００２９】
また、互いに隣接するノード点Ｎｉとノード点Ｎｊ間の流量は、次の式（６）で与えられる。
【００３０】
【数１】
2017041124.tif 000003
【００３１】
ここで、Ｃｉｊは流量係数、Ｚｉｊは弁開度（０〜１）である。
【００３２】
第１の流路の式（１）から式（３）を次の式（７）のように、ベクトルで表現する。
【００３３】
【数２】
2017041124.tif 000004
【００３４】
同様に、第２の流路の式（４）および式（５）を次の式（８）のように、ベクトルで表現する。
【００３５】
【数３】
2017041124.tif 000005
【００３６】
ただし、０と表示しているのは、同一要素数のゼロベクトルである。
【００３７】
また、図４に示す結合後の流路構成の場合は、次の式（９）のようなベクトル表示となる。
【００３８】
【数４】
2017041124.tif 000006
【００３９】
以上のように、流量バランスベクトル設定部１１は、第１の流路について式（７）、第２の流路について式（８）の形に設定する。なお、流量バランスベクトルのそれぞれの要素は、式（６）に示すように、各ノード点の圧力の関数である。
【００４０】
流量バランスベクトル設定部１１は、乗算部１１ａを有する。流量バランスベクトルの各要素を構成する各項の値のオーダーが異なる場合、後述する収束演算に支障をきたす可能性がある。このため、乗算部１１ａは、これらの値のオーダーの違いを緩和する演算が可能に構成されている。
【００４１】
具体的には、前述の式（９）に代えて、次の式（１０）に示すように、各要素に調整係数を乗じたものとする。
【００４２】
【数５】
2017041124.tif 000007
【００４３】
ここで、調整係数ｆ１ないしｆ５は、たとえば、それぞれ、次の式（１１）ないし（１５）のように流量係数の最大値の逆数とする。
【００４４】
【数６】
2017041124.tif 000008
【００４５】
【数７】
2017041124.tif 000009
【００４６】
【数８】
2017041124.tif 000010
【００４７】
【数９】
2017041124.tif 000011
【００４８】
【数１０】
2017041124.tif 000012
【００４９】
流量バランスベクトル記憶部２２は、流量バランスベクトル設定部１１で設定された流量バランスベクトルを記憶する。
【００５０】
変分比演算部１２は、流量バランスベクトル記憶部２２に記憶された流量バランスベクトルを用いて、変分比を演算する。
【００５１】
一般的に、関数Ｆ（ｘ）のヤコビアンＡ（ｘ）は、次の式（１６）により与えられる。
【００５２】
【数１１】
2017041124.tif 000013
【００５３】
したがって、σが十分小さいとすれば、次の式（１７）により近似できる。
【００５４】
【数１２】
2017041124.tif 000014
【００５５】
ここで、式（１７）の左辺のＡ（ｘ）ｐを変分比Ａｐと呼ぶこととする。変分比演算部１２は、流量バランスベクトル記憶部２２に記憶された流量バランスベクトルＦ（ｐ）（ただし、ｐは、後述する残差ベクトル）を用いて、変分比Ａｐを演算する。変分比Ａｐは、ヤコビアン行列Ａと残差ベクトルｐとの積に対応するものであるが、実際は、式（１７）の右辺に示す演算であり、ヤコビアン行列の演算は不要である。
【００５６】
なお、σの設定については、たとえば、次の式（１８）で求める。ここで、ｈは、おおよその計算誤差の２乗根である。また、ベクトルｘＴはベクトルｘの転置ベクトルである。
【００５７】
【数１３】
2017041124.tif 000015
【００５８】
変分比記憶部２３は、変分比演算部１２で算出された変分比Ａｐを記憶する。
【００５９】
収束演算部１３は、入力部が受け入れて流路系構成記憶部２１が記憶する各部のパラメータおよび境界条件に関するデータと、変分比記憶部２３に記憶された変分比Ａｐを用いて、各ノード点の圧力および各ノード点間の流量を算出する。具体的には、複数の圧力（圧力ベクトルＰ）の圧力推定値を設定して評価値を算出し、評価値に基づいて圧力推定値を順次更新する収束演算により、複数の圧力の値を算出し、その圧力から複数の流量の値を算出する。また収束演算部１３は、収束判定を行う。
【００６０】
ここで、収束計算は、たとえば、共役勾配法のアルゴリズムを用いることができる。いま、解くべき変数ベクトルをｘｋとし、一回の反復計算での修正量をｐｋ、ヤコビアンをＡ、反復計算での定数項をｂとする。ここで、ヤコビアン行列Ａは、ベクトルＦ（Ｐ）の各要素を圧力ベクトルＰの各要素により偏微分したものを要素とする行列である。
【００６１】
すなわち、次の関係式（１９）のようにおくと、式（２０）を解くこととなる。
【００６２】
【数１４】
2017041124.tif 000016
【００６３】
Ａｘｋ＝ｂ                          …（２０）
【００６４】
図６は、収束演算部における手順を示すフロー図である。繰り返し計算のアルゴリズムは、以下のようになる。
【００６５】
まず、最初は第１近似ｘ０を適当に設定する（ステップＳ１０１）。このｘ０を用いて、ｐ０＝ｒ０＝ｂ−Ａｘ０により残差ベクトルｐ０を求める（ステップＳ１０２）。
【００６６】
次に、ｑｋ＝Ａｐｋにより、ｑｋを求める（ステップＳ１０３）。次に、次の式（２１）により、αｋを算出する（ステップＳ１０４）。
【００６７】
【数１５】
2017041124.tif 000017
【００６８】
ここで、たとえば（ｐｋ，ｒｋ）は、ｐｋとｒｋとの内積を表す。
【００６９】
このαｋを用いて、次の式（２２）によりｘｋ＋１を算出する（ステップＳ１０５）。
ｘｋ＋１＝ｘｋ＋αｋｐｋ                        …（２２）
【００７０】
さらに、次の式（２３）によりｒｋ＋１を算出する（ステップＳ１０６）。
ｒＫ＋１
＝ｂ−Ａｘｋ＋１＝ｒｋ−αｋＡｐｋ＋１＝ｒｋ−αｋｑｋ＋１      …（２３）
【００７１】
また、次の式（２４）によりβｋを算出する（ステップＳ１０７）。
【００７２】
【数１６】
2017041124.tif 000018
【００７３】
次に、ｒｋ＋１および、このβｋを用いて、次の式（２５）により、新たな残差ベクトルｐｋ＋１を算出する（ステップＳ１０８）。
ｐｋ＋１＝ｒｋ＋１＋βｋｐｋ            …（２５）
【００７４】
次に、ステップＳ１０８で算出された残差ベクトルｐｋ＋１を前回得られた残差ベクトルｐｋと比較して収束判定を行う（ステップＳ１０９）。収束判定は、たとえば、両者の差のベクトルの絶対値｜ｐｋ＋１−ｐｋ｜の値を所定の値と比較することにより行う。収束していないと判定した場合（ステップＳ１０９  ＮＯ）は、ステップＳ１０３以下を繰り返す。収束したと判定した場合（ステップＳ１０９  ＹＥＳ）は、収束計算を終了する。
【００７５】
以上の収束計算においては、ヤコビアン行列Ａが単独に使用されることがなく、残差ベクトルとの積であるＡｐの形で使用される。具体的には、ステップＳ１０３、ステップＳ１０４、ステップＳ１０６、およびステップＳ１０７で用いられている。このＡｐは、前述の変分比演算部１２で述べた変分比Ａｐである。したがって、これらのステップでは、変分比演算部１２で変分比の演算が行われ、変分比記憶部２３に記憶された変分比Ａｐが、収束演算部１３で使用される。
【００７６】
また、以上の収束演算の各ステップでの演算結果は、収束状態記憶部２４に更新されながら記憶され、次の繰り返しステップで必要に応じて使用され、また、収束判定後は、バランスが成立した状態における各部の圧力、流量を与え、この結果を、出力部３２が出力する。
【００７７】
時間進展管理部１４は、時間的な応答を解析する際の時間ステップの管理を行う。端点条件が時間的に変化する場合に、流路系の各ノードの圧力等の状態量の時間的な変化を解析することになる。この場合、それぞれの時点において収束演算部１３による収束演算を行い、その時点を進めていくことになる。
【００７８】
この際、流路系の慣性が小さく、端点条件の変化に対して準静的に応答するのであれば、各時点の端点条件に対して収束演算をそれぞれ行うことのみで応答が算出できる。この際の時間進展管理部１４は、時間の進展とそれぞれの時点での端点条件の読み込ませの管理となる。
【００７９】
一方、流路系の慣性が無視できない場合は、過渡的な状態を記述した動特性方程式を解くことになる。この場合も、本実施形態により、演算可能である。すなわち、流路系の動特性方程式が式（２６）で与えられている場合、これを、式（２７）のように差分化して、流量を圧力の関数として式（２８）のように表すことができる。すなわち、式（２８）を前述の式（６）に代えて用いれば、すでに説明した手順と同様の手順で演算することができる。
【００８０】
【数１７】
2017041124.tif 000019
【００８１】
【数１８】
2017041124.tif 000020
【００８２】
【数１９】
2017041124.tif 000021
【００８３】
ここで、Ｉは慣性を表す係数、Δtはタイムステップ、Wijn+1は時刻n+1での流量を表す。
【００８４】
図７は、実施形態に係るシミュレーション方法の手順を示すフロー図である。まず、時間ｔをゼロに設定する（ステップＳ０１）。次に、入力部３１が、流路系の構成、各部のパラメータおよび境界条件に関するデータを入力として外部から受け入れ、流路系構成記憶部２１が流路構成その他を記憶する（ステップＳ０２）。次に、時間進展管理部１４が、時刻をΔｔだけ進める（ステップＳ０３）。次に、流量バランスベクトル設定部１１が流量バランスベクトルを設定し、流量バランスベクトル記憶部２２が記憶する（ステップＳ０４）。
【００８５】
次に、変分比演算部１２が変分比を演算し、収束演算部１３が収束演算を行う（ステップＳ１５）。収束演算部１３による収束計算の詳細は、図６に示した通りである。また、収束演算部１３は収束判定を行う（ステップＳ０６）。収束していないと判定された場合（ステップＳ０６  ＮＯ）は、ステップＳ０４以下を繰り返す。収束と判定された場合（ステップＳ０６  ＹＥＳ）は、時間進展管理部１４が時間の終了判定を行う（ステップＳ０７）。終了と判定されなかった場合（ステップＳ０７  ＮＯ）は、ステップＳ０３以下を繰り返す。終了と判定された場合（ステップＳ０７  ＹＥＳ）は、終了する。
【００８６】
図８は、実施形態に係るシミュレーション方法の結合後の系統についての計算の手順を示すフロー図である。これは、図２で示したような独立した流路が２つある場合の手順である。この場合は、図７のステップＳ０４ないしＳ０６のループを、図８に示すように、Ｓ０４、Ｓ１５、およびＳ１６のループの収束計算を行った後に、Ｓ２４ないしＳ２６のループの収束計算を行う。これを１つの収束演算で行うことは、ヤコビアン行列の次元が大きくなり、合理的ではない。
【００８７】
次に、本実施形態の効果を説明する。本実施形態の特徴を明確にするために、通常の方法と対比させながら説明する。
【００８８】
図９は、従来のシミュレーション方法の手順を示すフロー図である。本実施形態について示した図７との違いは、ステップＳ０５のヤコビアン行列を構成する部分である。
【００８９】
ニュートンラプソン法では、図２に示す体系についての圧力ベクトルｐの収束演算において、次のように、ヤコビアン行列の逆行列を求める必要がある。
【００９０】
【数２０】
2017041124.tif 000022
【００９１】
【数２１】
2017041124.tif 000023
【００９２】
ここで、前述の式（７）に示すように、ＦＡは第１の流路の流量バランスベクトル、ＰＡは第１の流路の圧力ベクトルであり、それぞれ３次元ベクトルである。したがって、ヤコビアン行列は、次の式（３１）に示す３×３行列となる。
【００９３】
【数２２】
2017041124.tif 000024
【００９４】
また、前述の式（８）に示すように、ＦＢは第２流路の流量バランスベクトル、ＰＢは第２流路の圧力ベクトルであり、それぞれ２次元ベクトルである。ヤコビアン行列は、次の式（３２）に示す２×２行列となる。
【００９５】
【数２３】
2017041124.tif 000025
【００９６】
ニュートンラプソン法では、前述のように、式（３１）および式（３２）の行列の逆行列を求める必要がある。
【００９７】
図４で示した第１の流路と第２の流路が結合している場合は、流路式のＦ（ｐ）は、前述の式（９）により与えられる５次元ベクトルである。この各要素を、５次元の圧力ベクトルの要素である端点を除くノード点の圧力、ｐ１、ｐ２、ｐ３、ｐ６、およびｐ７でそれぞれ偏微分すると、５×５行列のヤコビアンが得られる。ゼロ以外の要素は次の通りである。ただし、左辺の指標ｉｊは、ｉ行ｊ列の要素であることを表す。
【００９８】
【数２４】
2017041124.tif 000026
【００９９】
【数２５】
2017041124.tif 000027
【０１００】
【数２６】
2017041124.tif 000028
【０１０１】
【数２７】
2017041124.tif 000029
【０１０２】
【数２８】
2017041124.tif 000030
【０１０３】
【数２９】
2017041124.tif 000031
【０１０４】
【数３０】
2017041124.tif 000032
【０１０５】
【数３１】
2017041124.tif 000033
【０１０６】
【数３２】
2017041124.tif 000034
【０１０７】
【数３３】
2017041124.tif 000035
【０１０８】
【数３４】
2017041124.tif 000036
【０１０９】
【数３５】
2017041124.tif 000037
【０１１０】
【数３６】
2017041124.tif 000038
【０１１１】
以上のように、５×５行列と次元が大きくなり、収束計算の過程でその都度、５×５行列の逆行列を求める必要がある。一方、次元が大きくなる割には、５×５行列の２５の要素のうち、約半数の１２要素がゼロであり、比較的、疎な行列である。
【０１１２】
以上が、図９に示す従来の従来のシミュレーション方法の手順におけるニュートンラプソン法を、図２および図４の流路構成についてそれぞれ適用した具体例である。流路を結合した図４の流路構成の場合は、大きな次元のヤコビアン行列となり、その逆行列を収束過程でその都度演算する必要がある。
【０１１３】
一方、本実施形態においては、ヤコビアン行列そのものではなく、変分比Ａｐと呼んだＡ（ｘ）ｐを用いることにより、逆行列の演算を行う必要はなく、Ａ（ｘ）の次元の増大による影響が小さく、Ａ（ｘ）の次元の増大に容易に対応可能である。
【０１１４】
本実施形態においては、流路系構成記憶部２１が流量バランステーブルを有しており、流量構成の変更に対しても、流量バランステーブル内容が変更されることにより、流量バランスベクトル設定部１１により流量バランスベクトルが変更され、流路構成の変更が反映された結果が収束演算部１３での収束計算によって導かれる。
【０１１５】
また、流量バランスベクトルの各要素のそれぞれに調整係数が乗じられることにより、流量バランスベクトルの各要素の値が同程度のオーダーとなる。この結果、数値演算上、桁落ち等の問題を生ずることなく、収束演算部１３での収束計算が確実に行われる。
【０１１６】
以上のように、本実施形態によれば、プラントの流路網の計算について、プラント構成が一部変更された場合においても、シミュレーションにおける変更の手間を最小限とすることが可能となる。
【０１１７】
［その他の実施形態］
以上、本発明の実施形態を説明したが、実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。また、実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【０１１８】
１０…演算部、１１…流量バランスベクトル設定部、１１ａ…乗算部、１２…変分比演算部、１３…収束演算部、１４…時間進展管理部、２０…記憶部、２１…流路系構成記憶部、２２…流量バランスベクトル記憶部、２３…変分比記憶部、２４…収束状態記憶部、３１…入力部、３２…出力部、１００…シミュレーション装置
【図１】
2017041124.tif 000039
【図２】
2017041124.tif 000040
【図３】
2017041124.tif 000041
【図４】
2017041124.tif 000042
【図５】
2017041124.tif 000043
【図６】
2017041124.tif 000044
【図７】
2017041124.tif 000045
【図８】
2017041124.tif 000046
【図９】
2017041124.tif 000047