(57)【要約】【課題】効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる流体の数値計算方法及び装置及び制御プログラム及び記憶媒体を提供する。【解決手段】ＣＰＵ１１６により保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行う。【選択図】      図８
2005050063.tif 000002
【特許請求の範囲】
【請求項１】
保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うことを特徴とする流体の数値計算方法。
【請求項２】
保存形式の流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うことを特徴とする流体の数値計算方法。
【請求項３】
保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うように制御する制御手段を有することを特徴とする流体の数値計算装置。
【請求項４】
保存形式の流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うように制御する制御手段を有することを特徴とする流体の数値計算装置。
【請求項５】
請求項１または２記載の流体の数値計算方法をコンピュータにて実行するためのプログラムコードから成ることを特徴とするコンピュータ読み取り可能な制御プログラム。
【請求項６】
請求項５記載の制御プログラムを格納したことを特徴とする記憶媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、圧縮性を有する流体等の流体の数値計算方法及び装置及びその計算方法をコンピュータにて実行するための制御プログラム及びその制御プログラムを格納した記憶媒体に関する。
【０００２】
【従来の技術】
一般に、圧縮性を有する流体の数値解析手法は、大きくは、風上化を意識する解法と、格子系を工夫することにより風上化を意識しない解法とに分類することができる。
【０００３】
近年、効率的で安定性の高い圧縮性流体の風上化解法として、ＡＵＳＭ（Ａｄｖｅｃｔｉｏｎ  Ｕｐｓｔｒｅａｍ  Ｓｐｌｉｔｔｉｎｇ  Ｍｅｔｈｏｄ）スキームが提案されている（非特許文献１）。このＡＵＳＭスキームにおいては、非線形性の原因となる圧力を別扱いし、マッハ数に応じた圧力評価を行うために、状態方程式の関数形に応じた音速の評価を行っている。
【０００４】
また、風上化を意識しない解法としては、中心差分スキーム（セントラルスキーム）（非特許文献２）や、ＣＥ／ＳＥスキームが提案されている（非特許文献３）。中心差分スキームにおいては双対格子系を、ＣＥ／ＳＥスキームにおいては時空間における保存則を考慮するために、複数の格子系を採用し計算を行う。
【０００５】
【非特許文献１】
Ｍ．  Ｓ．  Ｌｉｕ  ａｎｄ  Ｃ．  Ｊ．  Ｓｔｅｆｆｅｎ，  Ｊｒ：Ａ  Ｎｅｗ  Ｆｌｕｘ  Ｓｐｌｉｔｔｉｎｇ  Ｓｃｈｅｍｅ，  Ｊｏｕｒｎａｌ  ｏｆ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ，  １０７，  ２３−３９（１９９３）
【非特許文献２】
Ｈ．  Ｎｅｓｓｙａｈｕ，  Ｅ．  Ｔａｄｍｏｒ，  “Ｎｏｎ−ｏｓｃｉｌｌａｔｏｒｙ  ｃｅｎｔｒａｌ  ｄｉｆｆｅｒｅｎｃｉｎｇ  ｆｏｒ  ｈｙｐｅｒｂｏｌｉｃ  ｃｏｎｓｅｒｖａｔｉｏｎ  ｌａｗｓ”  Ｊｏｕｒｎａｌ  ｏｆ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ，  ８７，  ４０８−４６３（１９９０）
【非特許文献３】
Ｓ．  Ｃ．  Ｃｈａｎｇ  “Ｔｈｅ  Ｍｅｔｈｏｄ  ｏｆ  Ｓｐａｃｅ−Ｔｉｍｅ  Ｃｏｎｓｅｒｖａｔｉｏｎ  Ｅｌｅｍｅｎｔ  ａｎｄ  Ｓｏｌｕｔｉｏｎ  Ｅｌｅｍｅｎｔ  −Ａ  Ｎｅｗ  Ａｐｐｒｏａｃｈ  ｆｏｒ  Ｓｏｌｖｉｎｇ  ｔｈｅ  Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ  ａｎｄ  Ｅｕｌｅｒ  Ｅｑｕａｔｉｏｎｓ”  Ｊｏｕｒｎａｌ  ｏｆ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ，  １１９，  ２９５−３２４（１９９５）
【０００６】
【発明が解決しようとする課題】
ＡＵＳＭスキームは、状態方程式の関数形に応じた音速の評価を行う必要があり、中心差分スキームやＣＥ／ＳＥスキームは、計算負荷が大きく、数値拡散が生じ易い複数の格子系を採用する必要がある。
【０００７】
本発明は、このような従来技術の有する事情に鑑みてなされたもので、その目的は、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる流体の数値計算方法及び装置及び制御プログラム及び記憶媒体を提供することである。
【０００８】
【課題を解決するための手段】
上記目的を達成するために本発明の流体の数値計算方法は、保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うことを特徴とする。
【０００９】
また、上記目的を達成するために本発明の流体の数値計算方法は、保存形式の圧縮性流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うことを特徴とする。
【００１０】
また、上記目的を達成するために本発明の流体の数値計算装置は、保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うように制御する制御手段を有することを特徴とする。
【００１１】
また、上記目的を達成するために本発明の流体の数値計算装置は、保存形式の圧縮性流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うように制御する制御手段を有することを特徴とする。
【００１２】
また、上記目的を達成するために本発明の制御プログラムは、前記流体の数値計算方法をコンピュータにて実行するためのプログラムコードから成ることを特徴とする。
【００１３】
更に、上記目的を達成するために本発明の記憶媒体は、前記制御プログラムを格納したことを特徴とする。
【００１４】
【発明の実施の形態】
以下、本発明の各実施の形態を図面に基づき説明する。
【００１５】
（第１の実施形態）
まず、本発明の第１の実施形態を図１乃至図６に基づき説明する。
【００１６】
下記式▲１▼は、一次元の場合の圧縮性を有する流体の基本方程式であり、Ｕは保存量の密度・運動量密度・単位体積当たりの全エネルギー、Ｆは流束である。ここで流束とは、下記式▲３▼で表される単位時間・単位体積当たりの量である。この方程式▲１▼は、保存量の時間変化と流束の収支の和が保存することを表している。
【００１７】
また、下記式▲２▼は、保存量の密度・運動量密度・単位体積当たりの全エネルギーＵを求めるための計算式である。
【００１８】
更に、下記式▲３▼は、流束Ｆを求めるための計算式である。
2005050063.tif 000003
【００１９】
図１は、本実施形態の特徴を示す図である。ここでは中心差分格子を採用しているが、流束を評価する際に、上記式▲３▼のように移流部分と圧力依存部分と散逸部分とに分割し、移流部分に関しては評価点における流速の方向に応じた風上差分を適用し、圧力依存部分に関しては中心差分によって計算した値を用いる。また、散逸部分については陽解法若しくは陰解法によって評価する。
【００２０】
図２は、一次元の場合の通常格子と中心差分格子とを示す図であり、同図においてｆｌｕｘはフラックス（流束）である。
【００２１】
図２の通常格子においてはセル中心とセル境界の位置は固定しているのに対し、中心差分格子においては双対格子系であり、セル中心とセル境界の位置が、１時間ステップ毎に半メッシュ移動するようになっている。また、通常格子においてはセル中心に保存量を定義し、セル境界において流束を評価する。また、中心差分格子においては保存量の定義位置と流束の評価位置とが同じである。
【００２２】
以下、本実施形態に係る圧縮性流体の数値計算方法による計算手順について、図３及び図４に基づき説明する。
【００２３】
図３及び図４は、本実施形態に係る圧縮性流体の数値計算方法による計算手順を示すフローチャートである。
【００２４】
まず、ステップＳ３０１において、保存量Ｕの初期設定を時刻ｎ及び時刻ｎ−１について行う（図４のステップＳ４０１）。
【００２５】
次に、ステップＳ３０２において、時刻ｎの内部エネルギー密度εを計算する（図４のステップＳ４０２）。
【００２６】
次に、ステップＳ３０３において、時刻ｎの圧力ｐを計算する（図４のステップＳ４０３）。
【００２７】
次に、ステップＳ３０４において、時刻ｎにおける流束Ｋと時刻ｎ＋１における流束の移流部分を計算する（図４のステップＳ４０４）。
【００２８】
移流部分の計算のための風上差分には、一次の風上差分だけでなく、高解像度スキーム（ＴＶＤ（Ｔｏｔａｌ  Ｖａｒｉａｔｉｏｎ  Ｄｉｍｉｎｉｓｈｉｎｇ）、ＥＮＯ（Ｅｓｓｅｎｔｉａｌｌｙ  Ｎｏｎ  Ｏｓｃｉｌｌａｔｏｒｙ）、ＰＰ（Ｐｏｓｉｔｉｖｅｌｙ  Ｐｒｅｓｅｒｖｉｎｇ）スキーム等）を用いることができる。
【００２９】
次に、ステップＳ３０５（図４のステップＳ４０５）において、保存量Ｕの更新を行い、時刻ｎ＋１の値を求める。ここで、流束の圧力依存部分に関しては中心差分を、移流部分については風上差分を、散逸部分に関しては時刻ｎ−１（陽解法）の値か、時刻ｎ＋１（陰解法）の値のいずれかを用いる。
【００３０】
本手法は、流束の圧力依存部分に関しては中心差分を適用するため、ＡＵＳＭスキームのように圧力評価のための音速を計算する必要がなく、状態方程式によらない汎用的な計算手法として利用することができる。
【００３１】
また、計算の安定条件が厳しく、数値振動が生じるような場合には、人工粘性項等を付与させることにより、計算の安定条件を満たすことが可能である。
【００３２】
図８は、本実施形態に係る圧縮性流体の数値計算方法を実行する計算機（流体の数値計算装置）の構成を示すブロック図である。
【００３３】
図８に示すように、この計算機は、複数のＣＰＵ１１６やメモリ１１８等からなる演算処理部、ハードディスク１１７やフロッピー（登録商標）ディスク１１０からなる記憶装置部、キーボード１１５やマウス１１４等からなる入力部、そしてディスプレイ１１２等からなる出力表示部等で構成されているシステムである。
【００３４】
この計算機は、図３及び４に示したフローチャートに対応するプログラムをメモリ１１８上にロードし、必要な計算領域を確保し、適当な方法で入力された領域形状や格子点座標及び物理量に基づき所定の演算処理を行い、その結果得られた領域上での物理量の値をハードディスク１１７等に書き込んで記憶させたり、ディスプレイ１１２上に表示させる。
【００３５】
図５は、圧縮性流体の代表的なベンチマークテストである衝撃波管問題を、本手法と中心差分スキームの一つであるＣＷＥＮＯスキームで計算した結果の一例を示す図である。
【００３６】
図６は、図５の一部を拡大した図である。
【００３７】
図５及び図６において、実線は解析解を、黒丸実線は本手法の計算結果を、三角点線はＣＷＥＮＯスキームの計算結果を示している。
【００３８】
図５及び図６にて明確なように、本手法の方が数値拡散が少なく、高い解像度が得られていることが分る。
【００３９】
以上詳述したように本実施の形態によれば、風上化解法の利点と非風上化解法の利点を採用し、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる。
【００４０】
また、本実施の形態によれば、計算スキームの変更無しに、理想気体の状態方程式のみならず、ｖａｎ  ｄｅｒ  Ｗａａｌｓ状態方程式等を扱うことができる。
【００４１】
（第２の実施形態）
次に、本発明の第２の実施形態を図７に基づき説明する。
【００４２】
下記式▲４▼は、二次元の場合の圧縮性流体の基本方程式であり、Ｕは保存量の密度・ｘ方向運動量密度・ｙ方向運動量密度・単位体積当たりの全エネルギー、Ｆはｘ方向流束、Ｇはｙ方向流束である。ここでｘ方向流束とは、下記式で表される単位時間・単位体積当たりの量である。この方程式▲４▼は、保存量の時間変化と流束の収支の和が保存することを表している。
【００４３】
また、下記式▲５▼は、保存量の密度・運動量密度・単位体積当たりの全エネルギーＵを求めるための計算式である。
【００４４】
また、下記式▲６▼は、ｘ方向流束Ｆを求めるための計算式である。
【００４５】
更に、下記式▲７▼は、ｙ方向流束Ｇを求めるための計算式である。
2005050063.tif 000004
【００４６】
図７は、二次元の場合の通常格子と中心差分格子の一例を示す図であり、同図（ａ）は通常格子を、同図（ｂ）は中心差分格子を示す。
【００４７】
この二次元は、上述した第１の実施形態における一次元の場合と同様に、通常格子においてはセル中心とセル境界の位置は固定しているのに対し、中心差分格子ではセル中心とセル境界の位置が、１時間ステップ毎に半メッシュ（ｘ方向だけでなく、ｙ方向にも）移動するようになっている。
【００４８】
また、この二次元は、一次元の場合と同様に、通常格子においてはセル中心に保存量を定義し、セル境界において流束を評価するが、中心差分格子においては保存量の定義位置と流束の評価位置とが同じである。
【００４９】
本発明を二次元の場合に適用するには、流束を評価する際に、上記式▲６▼、▲７▼のようにｘ方向流束とｙ方向流束とを移流部分と圧力依存部分と散逸部分とに分割し、移流部分に関しては評価点における流速の方向に応じた風上差分を適用し、圧力依存部分に関しては中心差分によって計算した値を用いれば良い。また、散逸部分については陽解法若しくは陰解法によって評価する。
【００５０】
これらの手続きは、ｙ方向の成分が付け加わったことを除けば、上述した第１の実施形態における一次元の場合と同様である。また、流体方程式系に、重力等の外力が付け加わった場合や、三次元の場合にも本手法は利用できる。
【００５１】
【発明の効果】
以上説明したように本発明によれば、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができるという効果を奏する。
【図面の簡単な説明】
【図１】本発明の第１の実施形態に係る流体の数値計算方法を示す概念図である。
【図２】本発明の第１の実施形態に係る流体の数値計算方法における通常格子とセントラル格子の概念図である。
【図３】本発明の第１の実施形態に係る流体の数値計算方法における計算手順を示すフローチャートである。
【図４】本発明の第１の実施形態に係る流体の数値計算方法における計算手順を示すフローチャートである。
【図５】本発明の第１の実施形態に係る流体の数値計算方法による計算結果と中心差分スキームによる計算結果の比較例を示す図である。
【図６】図５の一部を拡大した図である。
【図７】本発明の第２の実施形態に係る流体の数値計算方法における二次元の通常格子と中心差分格子の概念図である。
【図８】本実施形態に係る圧縮性流体の数値計算方法を実行する計算機（流体の数値計算装置）の構成を示すブロック図である。
【符号の説明】
１１０  フロッピー（登録商標）ディスク
１１２  ディスプレイ
１１４  マウス
１１５  キーボード
１１６  ＣＰＵ
１１７  ハードディスク
１１８  メモリ
【図１】
2005050063.tif 000005
【図２】
2005050063.tif 000006
【図３】
2005050063.tif 000007
【図４】
2005050063.tif 000008
【図５】
2005050063.tif 000009
【図６】
2005050063.tif 000010
【図７】
2005050063.tif 000011
【図８】
2005050063.tif 000012