(57)【要約】【課題】流体解析において，解析の精度を下げることなく速度場の数値振動を低減する計算方法および流体解析装置を提供する。【解決手段】運動量保存式による対流項の計算（Ｓ４）、拡散項の計算（Ｓ５）、圧力勾配項の計算（Ｓ６）に基づいて得られた速度分布に対し，人工粘性項の計算（Ｓ８）を行って速度の最新値を設定し（Ｓ９）、この値が質量保存則を満たすかを判定する。満たしたときは次の解析時刻の計算に移る。満たさないときは速度・圧力の修正計算を行って、質量保存則を満たすかを再判定する。人工粘性項の計算（Ｓ８）は２ｍ次微分オペレータ操作を行い，その結果に対し，係数の符号を変えて更に２ｍ次微分オペレータ操作を行う。【選択図】図２
2009169733.tif 000002
【特許請求の範囲】
【請求項１】
流体の速度を計算する流体解析装置であって，
前記流体解析装置は，一つもしくは複数の記憶装置と，一つもしくは複数の演算処理装置と，および前記記憶装置と前記演算装置との間でデータを通信するためのデータ通信手段を有しており，
前記記憶装置には，流体の速度を計算するための流体解析手順が格納されており，
前記速度は，前記流体解析手順に従い前記演算処理装置で計算され，前記記憶装置に格納されるようになっており，
前記流体解析手順は，流体の運動量保存式を計算する運動量保存式計算手順と人工粘性項計算手順を含んでおり，
前記運動量保存式計算手順は，流体の運動量保存式の対流項および圧力勾配項を計算する計算手順を含むもので，
前記人工粘性項計算手順は第１演算処理手順および第２演算処理手順からなり，前記流体解析手順の計算途中での流体の速度を暫定速度と呼ぶとき，
前記第１演算処理手順は，前記記憶装置に格納しておいた前記暫定速度を前記演算処理装置に呼び出し，前記演算処理装置において，ｍを自然数として，前記暫定速度の２ｍ次微分項に実係数Ｃを乗じた値に基づいて中間速度を計算し，計算された中間速度を前記記憶装置に格納する処理手順であり，
前記第２演算処理手順は，前記記憶装置に格納しておいた前記中間速度を前記演算処理装置に呼び出し，前記演算処理装置において，前記中間速度の２ｍ次微分項に実係数Ｃを乗じた値に基づいて人工粘性付与速度を計算し，計算された人工粘性付与速度を前記記憶装置に格納する処理手順であり，
前記人工粘性付与速度を用いてさらに前記流体解析手順に則り計算することで前記流体の速度を計算することを特徴とする流体解析装置。
【請求項２】
請求項１に記載の流体解析装置であって，
前記第１演算処理手順は，前記暫定速度の２ｍ次微分項に実係数Ｃを乗じた値を前記暫定速度に加算して中間速度を計算し，前記第２演算処理手順は，前記中間速度の２ｍ次微分項に実係数Ｃを乗じた値を前記中間速度から減算して人工粘性付与速度を計算する手順となる流体解析装置。
【請求項３】
請求項１に記載の流体解析装置であって，
前記第１演算処理手順は，前記暫定速度の２ｍ次微分項に実係数Ｃを乗じた値を前記暫定速度から減算して中間速度を計算し，前記第２演算処理手順は，前記中間速度の２ｍ次微分項に実係数Ｃを乗じた値を前記中間速度に加算して人工粘性付与速度を計算する手順となる流体解析装置。
【請求項４】
請求項１に記載の流体解析装置であって，
前記実係数Ｃの絶対値は，解析メッシュ幅を２ｍ乗した値を，円周率を２ｍ乗した値で除算した値以下で，かつ，０よりも大きい値になるように設定して，前記人工粘性項計算手順を実施することを特徴とした流体解析装置。
【請求項５】
請求項１に記載の流体解析装置であって，
前記自然数ｍの値を２以上に設定して前記人工粘性項計算手順を実施することを特徴とした流体解析装置。
【請求項６】
演算処理装置と記憶装置を有して流体の速度を計算する流体解析方法であって，
前記記憶装置に保存され、前記速度を計算するための流体解析計算手順は人工粘性項計算手順を含んでおり，
前記人工粘性項計算手順は第１演算処理手順および第２演算処理手順からなり，前記流体解析手順の計算途中での流体の速度を暫定速度と呼ぶとき，
前記第１演算処理手順は、ｍを自然数として，前記暫定速度の２ｍ次微分項に実係数Cを乗じて前記暫定速度に加算して中間速度を計算し，
前記第２演算処理手順は，前記中間速度の２ｍ次微分項に実係数Cを乗じて前記中間速度から減算して人工粘性付与速度を計算し，
計算された人工粘性付与速度を用い，さらに前記流体解析手順に則り計算することで流体の速度を計算することを特徴とする流体解析方法。
【請求項７】
演算処理装置と記憶装置を有して流体の速度を計算する流体解析方法であって，
前記記憶装置に保存され、前記速度を計算するための流体解析計算手順は人工粘性項計算手順を含んでおり，
前記人工粘性項計算手順は第１演算処理手順および第２演算処理手順からなり，前記流体解析手順の計算途中での流体の速度を暫定速度と呼ぶとき，
前記第１演算処理手順は，ｍを自然数として，前記暫定速度の２ｍ次微分項に実係数Cを乗じて前記暫定速度から減算して中間速度を計算し，
前記第２演算処理手順は，前記中間速度の２ｍ次微分項に実係数Ｃを乗じて前記中間速度に加算して人工粘性付与速度を計算し，
計算された人工粘性付与速度を用いて，さらに前記流体解析手順に則り計算することで流体の速度を計算することを特徴とする流体解析方法。
【請求項８】
流体解析の手順を記憶する記録媒体であって，
請求項６又は７に記載の流体解析方法の手順を読み出し可能に記憶したことを特徴とする記憶媒体。
【請求項９】
請求項１に記載の流体解析装置であって、計算された前記流体速度を出力手段から出力することを特徴とする流体解析装置。
【請求項１０】
請求項６に記載の流体解析方法であって、計算された前記流体速度を出力手段から出力することを特徴とする流体解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は流れ場などの速度を数値解析するための計算方法に関するものであり，特に数値振動の低減方法に関するものである。
【背景技術】
【０００２】
流体機器の設計や開発では，流体解析が利用されることが多い。流体解析手法の多くは，解析の対象となる領域を，空間的に有限個の検査体積，あるいは有限個の節点に分割して扱う。これを空間の離散化といい，一般的にはメッシュ分割，要素分割，セル分割などとよぶ。
【０００３】
空間の離散化に伴い，圧力や速度などの物理量も，検査体積あるいは節点に対応させて離散的に扱われる。このとき，物理量が定義された位置を物理量定義点とよぶ。
【０００４】
また質量の保存則，運動量の保存則などの物理法則に基づく基礎式は，有限個に離散化された物理量を用いて近似される。この近似式を離散化式とよぶ。速度や圧力などの物理量の空間分布は，離散化式を用いて計算する。
【０００５】
流体解析で得られた物理量の空間分布は，非物理的な振動分布となることがある。このような現象は数値振動，あるいは数値不安定とよばれている。数値振動した物理量の分布は計算誤差を含んでいる。また数値振動が生じた状態で計算を進めると，数値振動の振幅がさらに増大して計算誤差がさらに拡大することがある。流体解析の目的は，流れ場の物理量分布を正確に予測することなので，数値振動は好ましくない。
【０００６】
速度に対する数値振動の低減方法としては，運動量保存則の対流項の離散化に風上差分法を用いる方法がある。風上差分法には様々な種類の方法が提案されており，例えば非特許文献１に記載されている１次精度風上差分法や，非特許文献２に記載されているQUICK法などがある。
【０００７】
一般に対流項の離散化では，隣接しあう物理量定義点間の物理量分布を補間計算で近似し，その値を用いて物理量の対流項を計算する。風上差分法は，物理量の補間計算の際に，物理量定義点間の流れの方向を計算し，補間に用いる重み係数を，風上側の影響が大きくなるように設定する離散化方法である。なお，対流項の離散化式には，補間計算に伴う近似誤差が含まれる。すなわち対流項の離散化式は，本来の基礎式の対流項に，近似誤差に相当する項を加えたものになっている。近似誤差に相当する項は，打ち切り誤差項とよばれる。
【０００８】
風上差分法は，偶数次の微分項を打ち切り誤差項として持つ。一次精度風上差分法では2次微分項，QUICK法では4次微分項が打ち切り誤差項となっている。偶数次の微分項は，物理量の空間分布を滑らかにし，物理量の空間勾配を小さくする作用をもつ。2次微分項である運動量保存則の粘性項はその一例である。粘性項との作用の類似性から，風上差分法の偶数次微分項からなる打ち切り誤差項は，人工粘性項ともよばれる。また，人工粘性項により物理量の空間分布が滑らかになる作用は，数値拡散とよばれる。
【０００９】
風上差分法が数値振動の低減に有効なのは，人工粘性項の作用で，物理量の空間的な振動分布の振幅が小さくなるためである。
【００１０】
【非特許文献１】藤井孝蔵著，「流体力学の数値計算法」，東京大学出版会 (1994)
【非特許文献２】Leonard, B. P., "A stable and Accurate Convective Modeling Procedure Based on Quadratic Upstream Interpolation", Computer Methocs inApplied Mechanics and Engineering, vol.19, pp.59-98 (1979)
【非特許文献３】谷口伸行，「LES計算における数値不安定とその対策」，第14回数値流体力学シンポジウム B08-1 (2000)
【発明の開示】
【発明が解決しようとする課題】
【００１１】
しかし，風上差分法の人工粘性項はあくまでも誤差項である。そのため風上差分法を適用することで数値振動は低減するが，誤差の大きい結果しかえられないことが問題となっている。
【００１２】
一般に人工粘性項の微分の次数が低いほど誤差が大きくなる傾向がある。そのため人工粘性項が４次微分項であるQUICK法は，人工粘性項が２次微分項である1次精度風上差分法に比べて誤差が小さい。
【００１３】
しかし，非特許文献３で指摘されているように，高い解析精度が要求されるラージ・エディ・シミュレーション（ＬＥＳ）などの乱流場の流体解析では，ＱＵＩＣＫ法でも誤差が大きくて正しい結果が得られない。
【００１４】
そのためＬＥＳの計算では，運動量保存則の対流項の離散化に，２次精度中心差分法や４次精度中心差分法などの中心差分法が用いられることが多い。中心差分法の打ち切り誤差項は，数値拡散作用をもたない奇数次の微分項である。
【００１５】
しかし２次精度中心差分法や４次精度中心差分法などの中心差分法を用いる場合は，数値振動による誤差が大きくなるという問題があった。
【００１６】
本発明の目的は、上記従来技術の問題点に鑑み，数値振動を低減し，かつ解析誤差が小さい流体解析装置及び方法を提供することである。
【課題を解決するための手段】
【００１７】
上記目的を達成するための本発明は，流体の速度を計算する流体解析装置であって，前記速度を計算するための流体解析計算手順は，前記流体解析装置の計算機メモリもしくは前記流体解析装置とデータ通信が可能な記録媒体に保存されており，前記流体解析計算手順は人工粘性項計算手順を含んでおり，前記人工粘性項計算手順は第１演算処理手順および第２演算処理手順からなり，前記流体解析手順の計算途中での流体の速度を暫定速度と呼ぶとき，前記第１演算処理手順はｍを自然数として，前記暫定速度の２ｍ次微分項に実係数Cを乗じた値に基づいて中間速度を計算し，計算された中間速度を前記記憶装置に格納する処理手順であり，前記第２演算処理手順は，前記記憶装置に格納しておいた前記中間速度を前記演算処理装置に呼び出し，前記演算処理装置において，前記中間速度の２ｍ次微分項に実係数Ｃを乗じた値に基づいて人工粘性付与速度を計算し，計算された人工粘性付与速度を前記記憶装置に格納する処理手順であり，前記第２演算処理手順により計算された速度を用いてさらに前記流体解析手順に則り計算することで前記流体の速度を計算することを特徴とする。
【００１８】
前記第１演算処理手順は，前記暫定速度の２ｍ次微分項に実係数Ｃを乗じた値を前記暫定速度に加算して中間速度を計算し，前記第２演算処理手順は，前記中間速度の２ｍ次微分項に実係数Ｃを乗じた値を前記中間速度から減算して人工粘性付与速度を計算する手順となる。あるいは，前記第１演算処理手順は，前記暫定速度の２ｍ次微分項に実係数Ｃを乗じた値を前記暫定速度から減算して中間速度を計算し，前記第２演算処理手順は，前記中間速度の２ｍ次微分項に実係数Ｃを乗じた値を前記中間速度に加算して人工粘性付与速度を計算する手順となる。
【００１９】
また，前記実係数Ｃの絶対値を，解析メッシュ幅を２ｍ乗した値を，円周率を２ｍ乗した値で除算した値以下でかつ，０よりも大きい値になるように設定して前記人工粘性項計算手順を実施することを特徴とする。
【００２０】
また、前記自然数ｍの値を２以上に設定して前記人工粘性項計算手順を実施することを特徴とする。
【００２１】
本発明は，演算処理装置と記憶装置を有して流体の速度を計算する流体解析方法であって，前記記憶装置に保存され、前記速度を計算するための流体解析計算手順は人工粘性項計算手順を含んでおり，前記人工粘性項計算手順は第１演算処理手順および第２演算処理手順からなり，前記流体解析手順の計算途中での流体の速度を暫定速度と呼ぶとき，前記第１演算処理手順は、ｍを自然数として，前記暫定速度の２ｍ次微分項に実係数Cを乗じて前記暫定速度に加算（または減算）して中間速度を計算し，前記第２演算処理手順は，前記中間速度の２ｍ次微分項に実係数Cを乗じて前記中間速度から減算（または加算）して人工粘性付与速度を計算し，計算された人工粘性付与速度を用い，さらに前記流体解析手順に則り計算することで流体の速度を計算することを特徴とする。
【発明の効果】
【００２２】
本発明によれば，解析精度を下げることなく数値振動を低減する流体解析方式が提供され，高精度の流体解析の実施が可能となる。
【発明を実施するための最良の形態】
【００２３】
図１は，本発明の一実施形態であるところの流体解析装置の構成図である。流体解析装置１は，記憶装置としてハードディスク２およびメモリ３を備えており，演算処理装置としてＣＰＵ４を備えている。またハードディスク２とメモリ３とは，データ転送ライン５を介してデータ通信ができるようになっている。またメモリ３とＣＰＵ４は，データ転送ライン６を介してデータ通信ができるようになっている。ハードディスク２には，流体解析手順を収めたプログラムファイル７が格納されている。
【００２４】
流体解析装置１により流体の速度を計算する際は，プログラムファイル７に収められた流体解析手順に関するデータの一部もしくは全部が，データ転送ライン５を介してメモリ３に転送され，メモリ３に格納される。次いでメモリ３に格納された流体解析手順に従い，ＣＰＵ４において流体の速度が演算処理される。演算処理の際は，流体解析手順に従い随時，データ転送ライン６を介してメモリ３とＣＰＵ４との間でデータがやりとりされる。計算された流体の速度に関するデータはメモリ３に格納され，データ転送ライン５を介してハードディスク２に転送され，計算結果ファイルとしてハードディスク２に格納される。
【００２５】
図２に本発明による流体解析手順の一実施例を示す。本例は密度一定の非圧縮性流れを対象として流体の速度を非定常解析するものであり，乱流モデルにラージ・エディ・シミュレーション（ＬＥＳ）を採用している。基礎式は（１）式および（２）式である。
【００２６】
【数１】
2009169733.tif 000003
【００２７】
【数２】
2009169733.tif 000004
【００２８】
（１）式は運動量保存式であり，左辺第１項は時間変化項，右辺第1項は対流項，右辺第２項は圧力勾配項，右辺第３項は拡散項である。（２）式は質量保存式であり，左辺第1項は速度発散項である。
【００２９】
図２の流体解析手順はＳ１からＳ１４までの計算ステップから構成される。Ｓ１，Ｓ２，Ｓ３のステップは流体解析の準備段階のステップである。Ｓ４，Ｓ５，Ｓ６，Ｓ７，Ｓ８，Ｓ９，Ｓ１０，Ｓ１１，Ｓ１２，Ｓ１３，Ｓ１４のステップは非定常計算の時間反復ループを構成している。この流体解析手順では，時間刻みをΔｔとし，時間反復ループを１回おこなうたびに解析時刻をΔｔだけ進めるようになっている。
【００３０】
初めに計算条件入力ステップＳ１において,解析領域の形状データ、動粘性係数や密度などの物性値データ、時間刻みや収束判定値などの計算条件データの入力を行う。
【００３１】
図３に計算条件入力ステップＳ１でのデータフローを示す。ハードディスク２には形状データを収めた形状データファイル８，物性値データを収めた物性値データファイル９，計算条件データを収めた計算条件データファイル１０が格納されており，本ステップにおいてこれらのファイルに収められたデータをメモリ３に転送し，メモリ上に形状データ１１，物性値データ１２，計算条件データ１３として格納する。
【００３２】
次に初期値設定ステップＳ２において速度、圧力等の初期値を設定する。図４に初期値設定ステップＳ２でのデータフローを示す。ハードディスク２には初期値設定データを収めた初期値設定データファイル１４が格納されている。本ステップにおいて，初期値設定データファイル１４に収められた初期値設定データがメモリ３に転送され，初期値設定データ１５として格納される。次いでメモリ３から初期値設定データがＣＰＵ４に転送され，ＣＰＵ４において初期値設定処理演算が行われ，初期速度と初期圧力が計算される。これらのデータはＣＰＵ４からメモリ３に転送され，初期速度データ１６，初期圧力データ１７としてメモリ３に格納される。
【００３３】
次に解析時刻ゼロクリアステップＳ３において、解析時刻ｔの値を０に，タイムステップｎの値を０に設定する。
【００３４】
次に対流項計算ステップＳ４において、数１の対流項による加速度を計算する。この際，対流項の離散化には，２次精度中心差分法，あるいは４次精度中心差分法などの中心差分法を用いる。
【００３５】
図５に対流項計算ステップＳ４のデータフローを示す。メモリ３にはｎタイムステップでの速度データ１８が格納されている。タイムステップｎの値が０のときのｎタイムステップでの速度データ１８は，初期値設定ステップＳ２で計算された初期速度データ１６そのものである。
【００３６】
ｎタイムステップでの速度データ１８は，メモリ３からＣＰＵ４に転送され，ＣＰＵ４において（３）式に示す計算式を計算するのに用いられる。ＣＰＵ４で計算された対流項による加速度データはＣＰＵ４からメモリ３に転送され，対流項による加速度Ai,Cデータ１９としてメモリ３に格納される。
【００３７】
【数３】
2009169733.tif 000005
【００３８】
次に拡散項計算ステップＳ５において、（１）式の拡散項による加速度を計算する。図６に拡散項計算ステップＳ５のデータフローを示す。メモリ３に格納されているｎタイムステップでの速度データ１８，動粘性係数データ２０，サブグリッドスケール動粘性係数２１はＣＰＵ４に転送され，ＣＰＵ４において（４）式に示す計算式を計算するのに用いられる。ＣＰＵ４で計算された拡散項による加速度データはＣＰＵ４からメモリ３に転送され，拡散項による加速度Ai,Dデータ２２としてメモリ３に格納される。
【００３９】
【数４】
2009169733.tif 000006
【００４０】
次に圧力勾配項計算ステップＳ６において、（１）式の圧力勾配項による加速度を計算する。
【００４１】
図７に圧力勾配項計算ステップＳ６のデータフローを示す。メモリ３に格納されているｎタイムステップでの圧力データ２３，密度データ２４はＣＰＵ４に転送され，ＣＰＵ４において（５）式に示す計算式を計算するのに用いられる。ＣＰＵ４で計算された圧力勾配項による加速度データはＣＰＵ４からメモリ３に転送され，圧力勾配項による加速度Ai,Pデータ２５としてメモリ３に格納される。
【００４２】
【数５】
2009169733.tif 000007
【００４３】
次に速度更新ステップＳ７において，Ｓ４，Ｓ５，Ｓ６で計算した対流項、拡散項、圧力項による加速度を用いて速度を更新する。
【００４４】
図８に圧力勾配項計算ステップＳ６のデータフローを示す。メモリ３に格納されているｎタイムステップでの速度データ１８，対流項による加速度データ１９，拡散項による加速度データ２２，圧力勾配項による加速度データ２５はＣＰＵ４に転送され，ＣＰＵ４において（６）式に示す計算をするのに用いられる。ＣＰＵ４で計算された暫定速度データはＣＰＵ４からメモリ３に転送され，暫定速度データ２６としてメモリ３に格納される。
【００４５】
【数６】
2009169733.tif 000008
【００４６】
次に人工粘性項計算ステップＳ８において,人工粘性項を計算して速度を修正する。人工粘性項計算ステップＳ８は，図９に示すように３段階のステップからなり，初めに２ｍ次微分項の実係数Ｃを求める実係数計算ステップＳ８０，次に２ｍ次微分項計算（その１）ステップＳ８１，最後に２ｍ次微分項計算（その２）ステップＳ８２を行う。
【００４７】
図１０に実係数計算ステップＳ８０によるデータフローを示す。メモリ３には，メッシュ幅（Δｘ）データ２７，０より大きく１以下の定数ａに関する定数ａデータ２８，円周率データ２９，自然数ｍデータ３０が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（７）式を計算するのに用いられる。ＣＰＵ４で計算された実係数ＣのデータはＣＰＵ４からメモリ３に転送され，実係数Ｃデータ３１としてメモリ３に格納される。
【００４８】
【数７】
2009169733.tif 000009
【００４９】
次に２ｍ次微分項計算（その１）ステップＳ８１において，中間速度を計算する。図１１に２ｍ次微分項計算（その１）ステップＳ８１のデータフローを示す。メモリ３には暫定速度データ２６，実係数Ｃデータ３１，自然数ｍデータ３０が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（８）式を計算するのに用いられる。ＣＰＵ４で計算された中間速度のデータはＣＰＵ４からメモリ３に転送され，中間速度データ３２としてメモリ３に格納される。
【００５０】
【数８】
2009169733.tif 000010
【００５１】
次に２ｍ次微分項計算（その２）ステップＳ８２において，人工粘性付与速度を計算する。図１２に２ｍ次微分項計算（その２）ステップＳ８２のデータフローを示す。メモリ３には中間速度データ３２，実係数Ｃデータ３１，自然数ｍデータ３０が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（９）式を計算するのに用いられる。ＣＰＵ４で計算された人工粘性付与速度（ｕ♯）のデータはＣＰＵ４からメモリ３に転送され，人工粘性付与速度データ３３としてメモリ３に格納される。
【００５２】
【数９】
2009169733.tif 000011
【００５３】
次に最新値設定ステップ９において，タイムステップ（ｎ＋１）における速度および圧力を設定する。本ステップでは，（１０）式に従い，人工粘性付与速度（ｕ♯ｉ）を（ｎ＋１）タイムステップでの速度とし，また（１１）式の計算式に従い，ｎタイムステップの圧力を（ｎ＋１）タイムステップでの圧力としてメモリ３に格納する。
【００５４】
【数１０】
2009169733.tif 000012
【００５５】
【数１１】
2009169733.tif 000013
【００５６】
次に質量保存則満足度の判定ステップＳ１０において，（ｎ＋１）タイムステップでの速度が質量保存式を満足するか判定する。
【００５７】
図１３に質量保存則満足度の判定ステップＳ１０のデータフローを示す。メモリ３には（ｎ＋１）タイムステップでの速度データ３４，収束判定値データ３５が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（１２）式に示す不等式が満たされているか判定するのに用いられる。
【００５８】
【数１２】
2009169733.tif 000014
【００５９】
（１２）式が満たされていない場合には，次に速度・圧力修正計算ステップＳ１１の処理を行う。（１２）式が満たされている場合には，計算結果出力ステップＳ１２の処理を行う。
【００６０】
図１４に速度・圧力修正値計算ステップを示す。速度・圧力修正計算ステップＳ１１では，圧力修正値の計算ステップＳ１１１，速度修正計算ステップＳ１１２，圧力修正計算ステップ１１３を順次行う。
【００６１】
図１５に圧力修正値計算ステップＳ１１１のデータフローを示す。メモリ３には（ｎ＋１）タイムステップでの速度データ３４，時間刻みデータ３６，密度データ２４が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（１３）式を計算するのに用いられる。ＣＰＵ４で計算された圧力修正値ｐ'のデータはＣＰＵ４からメモリ３に転送され，圧力修正値３７としてメモリ３に格納される。
【００６２】
【数１３】
2009169733.tif 000015
【００６３】
図１６に速度修正計算ステップＳ１１２のデータフローを示す。メモリ３には（ｎ＋１）タイムステップでの速度データ３４，時間刻みデータ３６，密度データ２４，圧力修正値データ３７が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（１４）式を計算するのに用いられる。ＣＰＵ４では，（ｎ＋１）タイムステップでの速度を（１４）式の右辺の計算値で置換する。ＣＰＵ４では計算された（ｎ＋１）タイムステップでの速度データはＣＰＵ４からメモリ３に転送され，（ｎ＋１）タイムステップでの速度データ３４に上書きされメモリ３に格納される。
【００６４】
【数１４】
2009169733.tif 000016
【００６５】
図１７に圧力修正計算ステップＳ１１３のデータフローを示す。メモリ３には（ｎ＋１）タイムステップでの圧力データ３８，圧力修正値データ３７が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（１５）式を計算するのに用いられる。ＣＰＵ４では，（ｎ＋１）タイムステップでの圧力を（１５）式の右辺の計算値で置換する。ＣＰＵ４では計算された（ｎ＋１）タイムステップでの圧力データはＣＰＵ４からメモリ３に転送され，（ｎ＋１）タイムステップでの圧力データ３８に上書きされメモリ３に格納される。
【００６６】
【数１５】
2009169733.tif 000017
【００６７】
速度・圧力修正計算ステップＳ１１の実施後は，図２の質量保存則満足度の判定ステップＳ１０において，（ｎ＋１）タイムステップでの速度が質量保存式を満足するか再度判定する。判定ステップＳ１０において，（ｎ＋１）タイムステップでの速度が質量保存則を満たすと判断された場合は、計算結果出力ステップＳ１２において計算結果をモニタやハードディスク２などに出力する。
【００６８】
次に計算終了判定ステップＳ１３において，解析時刻ｔが解析終了時刻ｔｅｎｄを越えたかどうか判定する。解析時刻が解析終了時刻を越えたと判断されたときは、計算を終了する。解析終了時刻を越えていない場合は、解析時刻更新ステップＳ１４において速度データおよび圧力データの置き換えと，解析時刻およびタイムステップの更新を行う。
【００６９】
図１８に解析時刻更新ステップＳ１４の処理フローを示す。速度データ置換ステップＳ１４１，圧力データ置換ステップＳ１４２，解析時刻データ置換ステップＳ１４３，タイムステップデータ置換ステップＳ１４４からなる。
【００７０】
図１９に速度データ置換ステップＳ１４１のデータフローを示す。メモリ３にはｎタイムステップでの速度データ１８，（ｎ＋１）タイムステップでの速度データ３４が格納されており，このうち（ｎ＋１）タイムステップでの速度データ３４はＣＰＵ４に転送され，ＣＰＵ４において（１６）式を計算するのに用いられる。ＣＰＵ４で計算されたｎタイムステップでの速度データはＣＰＵ４からメモリ３に転送され，メモリ３上にすでに存在しているｎタイムステップでの速度データ１８に上書きされてメモリ３に格納される。
【００７１】
【数１６】
2009169733.tif 000018
【００７２】
図２０に圧力データ置換ステップＳ１４２のデータフローを示す。メモリ３にはｎタイムステップでの圧力データ２３，（ｎ＋１）タイムステップでの圧力データ３８が格納されており，このうち（ｎ＋１）タイムステップでの圧力データ３８はＣＰＵ４に転送され，ＣＰＵ４において（１７）式を計算するのに用いられる。ＣＰＵ４で計算されたｎタイムステップでの圧力データはＣＰＵ４からメモリ３に転送され，メモリ３上にすでに存在しているｎタイムステップでの圧力データ２３に上書きされてメモリ３に格納される。
【００７３】
【数１７】
2009169733.tif 000019
【００７４】
図２１に解析時刻データ置換ステップＳ１４３のデータフローを示す。メモリ３には解析時刻データ３９，時間刻みデータ３６が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（１８）式を計算するのに用いられる。ＣＰＵ４では解析時刻ｔが置換計算され，その結果はＣＰＵ４からメモリ３に転送され，メモリ３上にすでに存在している解析時刻データ３９に上書きされてメモリ３に格納される。
【００７５】
【数１８】
2009169733.tif 000020
【００７６】
図２２にタイムステップデータ置換ステップＳ１４４のデータフローを示す。メモリ３にはタイムステップｎに関するタイムステップデータ４０が格納されており，本データはＣＰＵ４に転送され，ＣＰＵ４において（１９）式を計算するのに用いられる。ＣＰＵ４ではタイムステップｎが置換計算され，その結果はＣＰＵ４からメモリ３に転送され，メモリ３上にすでに存在しているタイムステップデータ４０に上書きされてメモリ３に格納される。以上の解析時刻更新ステップＳ１４の実施後は，対流項計算ステップＳ４を実施する。
【００７７】
【数１９】
2009169733.tif 000021
【００７８】
以上に示した流体解析手順にしたがって流体解析を行えば，解析精度を下げることなく，速度の数値振動を低減することができる。それは，人工粘性項計算ステップＳ８が，（１）式の運動量保存式に，４ｍ次微分項からなる人工粘性項を付与することに相当するからである。ここではまず，人工粘性項計算ステップＳ８が数１の運動量保存式に，４ｍ次微分項からなる人工粘性項を付与することに相当する理由を説明する。
【００７９】
（２０）式に示す計算式は２ｍ次微分オペレータ（その１）の定義式であり，（８）式を任意のスカラーｆに対して一般化して表記したものである。また（２１）式に示す計算式は２ｍ次微分オペレータ（その２）の定義式であり，（９）式を任意のスカラーｆに対して一般化して表記したものである。さらに（２２）式に示す計算式は，スカラーｆに２ｍ次微分オペレータ（その１）を作用させ，その結果に２ｍ次微分オペレータ（その２）を作用させた結果を示す式である。（２２）式は，ｆの４ｍ次微分項を，ｆから減じる式となっていることがわかる。
【００８０】
【数２０】
2009169733.tif 000022
【００８１】
【数２１】
2009169733.tif 000023
【００８２】
【数２２】
2009169733.tif 000024
【００８３】
人工粘性項計算ステップＳ８の計算において，２ｍ次微分項計算（その１）ステップＳ８１を行った後，２ｍ次微分項の計算（その２）ステップＳ８２を行うことは，暫定速度に２ｍ次微分オペレータ（その１）を作用させ，その結果に２ｍ次微分オペレータ（その２）を作用させたことに等しい。そのため人工粘性項計算ステップＳ８での計算は，４ｍ次微分項からなる人工粘性項を計算したことに相当する。
【００８４】
つぎに実係数Ｃの値の決め方について説明する。ここで４ｍ次微分オペレータＨを（２３）式で定義しておく。４ｍ次微分オペレータＨは，２ｍ次微分オペレータ（その１）の結果に２ｍ次微分オペレータ（その２）を作用させるのと等価なオペレータである。（２４）式のように，スカラーｆに４ｍ次微分オペレータを作用させた結果をｆ*とする。
【００８５】
【数２３】
2009169733.tif 000025
【００８６】
【数２４】
2009169733.tif 000026
【００８７】
（２４）式を対象として，フォン・ノイマンの安定性解析手法を適用する。スカラーｆの空間分布を，（２５）式で与える。（２５）式は，数値振動した計算結果を模擬するものであり，振幅がｇ，波数がｋの振動分布となっている。（２５）式を（２４）式に代入して，（２６）式を得る。ｆ*は，振幅が（１−C2k4m）ｇ，波数がｋの振動分布となっていることがわかる。
【００８８】
【数２５】
2009169733.tif 000027
【００８９】
【数２６】
2009169733.tif 000028
【００９０】
（２４）式に関する振幅の増幅率Ｇは，（２７）式で与えられる。増幅率Ｇが１以上だと，（２４）式を繰り返し計算した際に数値振動の振幅が増大してしまう。そこで数値振動を低減するためには，増幅率Ｇが（２８）式に示す不等式を満たす必要がある。（２８）式の不等式に（２７）式を代入し，実係数Ｃに関する（２９）式に示す不等式を得る。
【００９１】
【数２７】
2009169733.tif 000029
【００９２】
【数２８】
2009169733.tif 000030
【００９３】
【数２９】
2009169733.tif 000031
【００９４】
流体解析で生じる数値振動の最小の波長はメッシュ幅Δｘの２倍なので，流体解析で計算可能な波数の範囲は，（３０）式に示す不等式の範囲となる。（３０）式の条件より，数値振動を低減するための条件は（３１）式に示す不等式で与えられる。以上のことから，実係数Ｃの値は，（３１）式の不等式を満たすように与える。本実施例においては，実係数計算ステップＳ８０における実係数Ｃの計算は，（３１）式を満たすようになっている。
【００９５】
【数３０】
2009169733.tif 000032
【００９６】
【数３１】
2009169733.tif 000033
【００９７】
次に，自然数ｍの決め方について説明する。図２３に，人工粘性項の微分の次数４ｍと増幅率Ｇの関係を示す。ここでは実係数Ｃを（３２）式のように与えた。ｍの値によらず，増幅率Ｇは，波数が０のときに１となっており，波数が最大波数のときに０となっている。
ただし０＜ｋ＜ｋmaxの範囲の増幅率Ｇの値は，ｍの値が大きいほど１に近い値となっている。
【００９８】
【数３２】
2009169733.tif 000034
【００９９】
一般に流体解析で問題となるのは,速度に関する最大波数の数値振動である。この数値振動の波長はメッシュ幅の２倍であり，そのときの速度分布はチェッカーボード状，あるいはジグザグ状となる。また最大波数よりも波数の小さい振動は，多くの場合は数値振動ではなく，物理的に正しい振動である。
【０１００】
そこで,人工粘性項に求められる特性は，最大波数の振動を選択的に減衰させ，それよりも波数の小さい振動は減衰させないことである。すなわち最大波数での増幅率は０に近く，またそれより低波数での増幅率は１に近いことが求められる。この特性は，ｍの値を大きくするほど満たされるようになるので，ｍの値はなるべく大きめに設定することが望ましい。
【０１０１】
ｍ＝１の場合の人工粘性項の微分の次数は４次で，誤差の大きいQUICK法の人工粘性項と同じ次数である。そのためｍ＝１とすることは好ましくない。そのためｍには，２以上の値を設定する。自然数ｍを２以上の値として設定して人工粘性項計算ステップＳ８を実施すると，８次以上の人工粘性項を計算したことになる。
【０１０２】
図２３より,８次以上の人工粘性項は，ＱＵＩＣＫ法の人工粘性項である４次の人工粘性項に比べて低波数域での減衰が少なく，高波数の振動を選択的に減衰させることがわかる。この特性により，本発明による流体解析手順では，物理的に意味のある低波数域の振動は減衰させず，高波数域に発生する数値振動を選択的に減衰することができるので，解析精度を下げることなく，速度の数値振動を低減すること可能となる。
【０１０３】
なお（３３）式に示すように，スカラーｆに２ｍ次微分オペレータ（その２）を先に作用させ，その結果に２ｍ次微分オペレータ（その１）を作用させた結果は，スカラーｆに２ｍ次微分オペレータ（その１）を先に作用させ，その結果に２ｍ次微分オペレータ（その２）を作用させた結果と等しい。そこで人工粘性項計算ステップＳ８の代わりに，図２４に示す人工粘性計算ステップＳ８'を用いてもかまわない。
【０１０４】
【数３３】
2009169733.tif 000035
【０１０５】
人工粘性項計算ステップＳ８'は３段階のステップからなり，初めに２ｍ次微分項の実係数Ｃを求める実係数計算ステップＳ８０，次に２ｍ次微分項計算（その２'）ステップＳ８２'，最後に２ｍ次微分項計算（その１'）ステップＳ８１'を行う。
【０１０６】
２ｍ次微分項計算（その２'）ステップＳ８２'では,（３４）式を用いて中間速度を計算する。図２５に２ｍ次微分項計算（その２'）ステップＳ８２'のデータフローを示す。メモリ３には暫定速度データ２６，実係数Ｃデータ３１，自然数ｍデータ３０が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（３４）式を計算するのに用いられる。ＣＰＵ４で計算された中間速度のデータはＣＰＵ４からメモリ３に転送され，中間速度データ３２としてメモリ３に格納される。
【０１０７】
【数３４】
2009169733.tif 000036
【０１０８】
２ｍ次微分項計算（その１'）ステップＳ８１'では，(３５)式を用いて人工粘性付与速度を計算する。図２６に２ｍ次微分項計算（その１'）ステップＳ８１'のデータフローを示す。メモリ３には中間速度データ３２，実係数Ｃデータ３１，自然数ｍデータ３０が格納されており，これらのデータはＣＰＵ４に転送され，ＣＰＵ４において（３５）式を計算するのに用いられる。ＣＰＵ４で計算された人工粘性付与速度（ｕ♯）のデータはＣＰＵ４からメモリ３に転送され，人工粘性付与速度データ３３としてメモリ３に格納される。
【０１０９】
【数３５】
2009169733.tif 000037
【０１１０】
またＳ８の人工粘性項の計算は，必ずしもＳ７の速度更新の次に行う必要はない。例えば図２７のフローチャートに示すように，Ｓ８はＳ４の対流項の計算の直前に行っても良い。あるいは図２８のフローチャートに示すように，Ｓ１０の質量保存則の判定後にＳ８を行っても良い。
【０１１１】
上述した流体解析方法は、流体解析方法のステップをコンピュータに実行させるプログラムとして記録媒体に記録することができる。プログラムをコンピュータにインストールすることで、流体解析方法のステップをコンピュータに実行させる流体解析装置とすることができる。
【０１１２】
上述した流体解析装置および流体解析方法は、ボイラー、ガスタービン、自動車の流体解析に用いることができる。流体解析の対象は、例えば、ボイラーの燃焼ガス、供給燃料、供給空気、生成蒸気、ガスタービンの燃焼ガス、供給燃料、供給空気、自動車の燃焼ガス、供給燃料、供給空気などである。
【０１１３】
上述した流体解析装置および流体解析方法は、計算された流体速度を出力手段から出力する。出力手段はコンピュータの外部出力インターフェイスなどである。計算された流体速度は、表示装置へ出力されたり、他のコンピュータへ出力されたり、流体解析装置内のハードディスクや他のハードディスクへ出力されたり、流体解析装置内のメモリや他のメモリへ出力される。出力された流体解析結果は、別の流体解析に使用したり、制御対象の操作量を算出する為に使用することができる。例えば、ボイラーの燃焼ガスの流体解析結果をメモリに出力し、メモリへ記録した後、生成蒸気量の流体解析を行ったり、ボイラーの供給空気の操作量を算出する為、メモリ内の燃焼ガスの流体解析結果を参照し、使用することができる。
【図面の簡単な説明】
【０１１４】
【図１】本発明の実施の一形態である流体解析装置の構成図。
【図２】本発明の実施の一形態である流体解析の手順を示すフローチャート。
【図３】入力ステップＳ１でのデータフロー図。
【図４】初期値設定ステップＳ２でのデータフロー図。
【図５】対流項計算ステップＳ４でのデータフロー図。
【図６】拡散項計算ステップＳ５でのデータフロー図。
【図７】圧力項計算ステップＳ６でのデータフロー図。
【図８】速度の更新ステップＳ７でのデータフロー図。
【図９】人工粘性項計算ステップＳ８の手順を示すフローチャート。
【図１０】実係数の計算ステップＳ８０でのデータフロー図。
【図１１】２ｍ次微分項の計算（その１）ステップＳ８１でのデータフロー図。
【図１２】２ｍ次微分項の計算（その２）ステップＳ８２でのデータフロー図。
【図１３】質量保存則を満たすかの判定ステップＳ１０でのデータフロー図。
【図１４】速度・圧力修正計算ステップＳ１１の手順を示すフローチャート。
【図１５】圧力修正値の計算ステップＳ１１１でのデータフロー図。
【図１６】速度修正計算ステップＳ１１２でのデータフロー図。
【図１７】圧力修正計算ステップＳ１１３でのデータフロー図。
【図１８】解析時刻更新ステップＳ１４の手順を示すフローチャート。
【図１９】速度データ置換ステップＳ１４１でのデータフロー図。
【図２０】圧力データ置換ステップＳ１４２でのデータフロー図。
【図２１】解析時刻データ置換ステップＳ１４３でのデータフロー図。
【図２２】タイムステップデータ置換ステップＳ１４４でのデータフロー図。
【図２３】人工粘性項の特性の説明図。
【図２４】人工粘性項計算ステップの変形例を示すフローチャート。
【図２５】２ｍ次微分項の計算（その２'）ステップＳ８２'でのデータフロー図。
【図２６】２ｍ次微分項の計算（その１'）ステップＳ８１'でのデータフロー図。
【図２７】他の実施形態による流体解析の手順を示すフローチャート。
【図２８】更に他の実施形態による流体解析の手順を示すフローチャート。
【符号の説明】
【０１１５】
１…流体解析装置、２…ハードディスク、３…メモリ、４…ＣＰＵ、５…データ転送ライン、６…データ転送ライン、７…プログラムファイル、８…形状データファイル、９…物性値データファイル、１０…計算条件データファイル、１１…形状データ、１２…物性値データ、１３…計算条件データ、１４…初期値設定データファイル、１５…初期値設定データ、１６…初期速度データ、１７…初期圧力データ、１８…ｎタイムステップでの速度データ、１９…対流項による加速度Ai,Cデータ、２０…動粘性係数データ、２１…サブグリッドスケール動粘性係数データ、２２…拡散項による加速度Ai,Dデータ、２３…ｎタイムステップでの圧力データ、２４…密度データ、２５…圧力勾配項による加速度Ai,Pデータ、２６…暫定速度データ、２７…メッシュ幅データ、２８…定数ａデータ、２９…円周率データ、３０…自然数ｍデータ、３１…実係数Ｃデータ、３２…中間速度データ、３３…人工粘性付与速度データ、３４…（ｎ＋１）タイムステップでの速度データ、３５…収束判定値データ、３６…時間刻みデータ、３７…圧力修正値データ、３８…（ｎ＋１）タイムステップでの圧力データ、３９…解析時刻データ、４０…タイムステップｎデータ。
【図１】
2009169733.tif 000038
【図２】
2009169733.tif 000039
【図３】
2009169733.tif 000040
【図４】
2009169733.tif 000041
【図５】
2009169733.tif 000042
【図６】
2009169733.tif 000043
【図７】
2009169733.tif 000044
【図８】
2009169733.tif 000045
【図９】
2009169733.tif 000046
【図１０】
2009169733.tif 000047
【図１１】
2009169733.tif 000048
【図１２】
2009169733.tif 000049
【図１３】
2009169733.tif 000050
【図１４】
2009169733.tif 000051
【図１５】
2009169733.tif 000052
【図１６】
2009169733.tif 000053
【図１７】
2009169733.tif 000054
【図１８】
2009169733.tif 000055
【図１９】
2009169733.tif 000056
【図２０】
2009169733.tif 000057
【図２１】
2009169733.tif 000058
【図２２】
2009169733.tif 000059
【図２３】
2009169733.tif 000060
【図２４】
2009169733.tif 000061
【図２５】
2009169733.tif 000062
【図２６】
2009169733.tif 000063
【図２７】
2009169733.tif 000064
【図２８】
2009169733.tif 000065