(57)【要約】  【課題】  非圧縮性流体の非定常挙動を記述した連立方程式の解が収束するまでの時間を飛躍的に減少する非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体を提供する。  【解決手段】  解析点での流速ベクトル及び圧力の初期値を受け付け、解析点での流速ベクトル近似解を算出し、流速ベクトル近似解の発散の絶対値の最大値が所定値より大きい場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出した圧力修正量と記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいか再度判断し、流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さい場合、流速ベクトル近似解及び圧力近似解を出力する。【選択図】        図２
2005240889.tif 000002
【特許請求の範囲】
【請求項１】
ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置において、
求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶する記憶手段と、
前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、
受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、
前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、
該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、
該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、
算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、
算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、
前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段と
を備えることを特徴とする非圧縮性流体の非定常解析装置。
【請求項２】
前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、
該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、
該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段と
を備えることを特徴とする請求項１に記載の非圧縮性流体の非定常解析装置。
【請求項３】
ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求めるコンピュータを用いた非圧縮性流体の非定常解析方法において、
前記コンピュータは、
求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶し、
前記解析点での流速ベクトル及び圧力の初期値を受け付け、
受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出し、
前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断し、
前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出し、
算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、
算出した新たな圧力修正量を算出回数と対応付けて記憶し、
算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を用いて流速ベクトル近似解を算出し、再度流速ベクトル近似解の発散を算出し、算出した発散の絶対値の最大値が所定値より小さいか否か判断し、
前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力することを特徴とする非圧縮性流体の非定常解析方法。
【請求項４】
前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出し、
算出した残差の二乗和の平方根を最小とするパラメータを算出し、
算出したパラメータを用いて、記憶してある圧力修正量を線形に結合することを特徴とする請求項３に記載の非圧縮性流体の非定常解析方法。
【請求項５】
ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムにおいて、
求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、
前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、
受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、
前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、
該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、
該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、
算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、
算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、
前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段と
を備えることを特徴とするコンピュータプログラム。
【請求項６】
前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、
該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、
該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段と
を備えることを特徴とする請求項５に記載のコンピュータプログラム。
【請求項７】
ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムを記録した記録媒体において、
求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、
前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、
受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、
前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、
該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、
該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、
算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、
算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、
前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段と
を備えることを特徴とする記録媒体。
【請求項８】
前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、
該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、
該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段と
を備えることを特徴とする請求項７に記載の記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非圧縮性流体の非定常挙動を有限要素法、差分法、有限体積法等を用いて有限な解析点に離散化して解析する非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体に関する。
【背景技術】
【０００２】
従来の非圧縮性流体の非定常挙動を解析する方法について、層流を例に挙げて説明する。なお、例えば乱流であっても、非定常挙動の解析方法として本質的な相違点は生じない。
【０００３】
非圧縮性流体の非定常挙動を解析する場合、解を求めるべき方程式は、（数１）、（数２）で表されるナビエ・ストークスの運動方程式及び流体の連続性条件式である。（数１）、（数２）において、ベクトルｕは流体の流速ベクトルを、ｐは圧力を流体密度で除算した値を、ｔは時間を、νは動粘性係数を、それぞれ示す。
【０００４】
【数１】
2005240889.tif 000003
【０００５】
【数２】
2005240889.tif 000004
【０００６】
（数１）、（数２）の解を求める手段としては、非特許文献１に開示されているＳＭＡＣ法、非特許文献２に開示されているＦｒａｃｔｉｏｎａｌＳｔｅｐ法等が良く知られている。しかし、いずれの方法においても、解の収束条件に流体の連続条件式が含まれておらず、流体の状況が急変した場合においては、解の精度が保証できないという問題点があった。
【０００７】
そこで、例えば非特許文献３に開示されているＳＯＬＡ法、非特許文献４に開示されているＴＡＫＥＭＩＴＳＵの方法が提案されている。図６は、非特許文献３又は非特許文献４に開示されている従来の解析方法に従って、コンピュータを用いて処理した場合の手順を示すフローチャートである。
【０００８】
図６では、まず流速ベクトルｕ及び圧力ｐの初期値が入力され（ステップＳ６０１）、流速ベクトルｕの初期近似解を（数３）に基づいて算出する（ステップＳ６０２）。なお（数３）において、ベクトルの上部に付された波線（〜）は近似解であることを示しており、ｎは所定の時刻における状態であることを示している。
【０００９】
【数３】
2005240889.tif 000005
【００１０】
（数１）から（数３）を差し引くことにより、（数４）を求める。（数４）では、右辺は所定の時刻から単位計算時間Δｔが経過した後の圧力から、所定の時刻での圧力を引いた圧力修正量δｐの勾配、すなわちグラーディエント（gradient）δｐを、左辺は流速ベクトルｕの単位計算時間Δｔ当たりの変化を、それぞれ示している。
【００１１】
【数４】
2005240889.tif 000006
【００１２】
次に、（数４）の両辺に対して発散を取り、（数２）を代入して得ることができる（数５）に示すように、圧力修正量δｐの勾配の発散、すなわち流速ベクトルｕの発散を単位時間で除した値を算出する。
【００１３】
【数５】
2005240889.tif 000007
【００１４】
ｕn+1及びｐn+1を求めるために、任意の回数ｋだけ反復演算する。すなわち、カウンタｋの初期値を０（ゼロ）とし（ステップＳ６０３）、流速ベクトルｕの初期値を（数３）で算出した流速ベクトルｕの近似解とし（ステップＳ６０４）、全ての解析点における流速ベクトルｕkについて発散divｕkの絶対値の最大値が閾値であるεより小さいか否かを判断する（ステップＳ６０５）。流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより大きいと判断した場合（ステップＳ６０５：ＮＯ）、（数５）より圧力修正量δｐを算出し（ステップＳ６０６）、緩和係数ωを用いて、（数６）のように所定の時刻における流速ベクトルｕk及び圧力ｐkを算出して（ステップＳ６０７）、所定の時刻における流速ベクトルｕkを（数３）で算出した流速ベクトルｕの近似解とし（ステップＳ６０８）、カウンタｋを１インクリメントして（ステップＳ６０９）、ステップＳ６０４へ戻る。
【００１５】
【数６】
2005240889.tif 000008
【００１６】
流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより小さいと判断した場合（ステップＳ６０５：ＹＥＳ）、（数７）に示すように、求めた流速ベクトルｕk、圧力ｐkを、単位計算時間Δｔ経過後の流速ベクトルｕn+1、圧力ｐn+1として算出する（ステップＳ６１０）、上述した処理を目的の時刻ｔmaxになるまで繰り返すことで所定の時刻における流体の流速ベクトルｕ及び圧力ｐを求めることができる。
【００１７】
【数７】
2005240889.tif 000009
【００１８】
本方法によれば、流体の連続の式を満足しているか否かを、各次官ステップ毎に判定しており（ステップＳ６０５）、全ての時間ステップにおいて、必ず連続の式を満足した解を求めることが可能となる。
【非特許文献１】エイ．エイ．アムスデン（A.A.Amsden），エフ．エイチ．ハーロー（F.H.Harlow）著、「ＳＭＡＣ法：非圧縮性流体の数値計算技術」（The SMAC method: A numerical technique for calculating incompressible fluid flow）、ロスアラモス科学研究所（Los Alamos Sci. Lab.）、LA-4370、１９７０年
【非特許文献２】河村洋、土方邦夫編、「熱と流れのシミュレーション」、丸善株式会社、１９９５年、ｐ．５
【非特許文献３】シー．ダブリュ．ハート（C.W.Hirt）,ビー．ディー．ニクルス（B.D.Nichls）エヌ．シー．ロメロ（N.C.Romero）著、「ＳＯＬＡ：瞬間流体の数値解析アルゴリズム」（SOLA: A Numerical Solution Algorithm for Transient Fluid Flows）、ロスアラモス科学研究所（Los Alamos Sci. Lab.）、LA-5852、１９７５年
【非特許文献４】エヌ．タケミツ（N.Takemitsu）著、「非圧縮性流体の差分解法」（Finite Difference Method to Solve Incompressible Fluid Flow）、ジャーナルオブコンピュテーショナルフィジックス（Journal of Computational Physics）、vol.61、１９８５年、ｐ．４９９
【発明の開示】
【発明が解決しようとする課題】
【００１９】
しかし、上述したような非圧縮性流体の非定常挙動を解析する手法では、非圧縮性流体の非定常過程を厳密に追跡することは可能であるが、各時間ステップ毎の流速ベクトルと圧力とを求める収束処理を高速に行うことは困難であり、所望の時刻ｔmaxにおける流速ベクトル及び圧力の解を求めるまでに相当の時間を要するという問題点が残されている。
【００２０】
各時間ステップ毎の流速ベクトルと圧力とを求める処理を高速に行うべく、連立方程式の高速解法である残差切除法の原理を利用することも考えられている。最も簡単な方法として、流速ベクトル修正量を求めるステップに残差切除法を適用することが考えられる。
【００２１】
しかし、流速ベクトル修正量は、反復解法を適当な回数だけ反復して求めた近似的な値で十分であり、残差切除法を用いて完全に解を求める場合、たとえ残差切除法を適当な回数で打ち切り近似解を求めるようにした場合であっても、逆に解析時間が増加することは明らかであり好ましくない。
【００２２】
本発明は斯かる事情に鑑みてなされたものであり、非圧縮性流体の非定常挙動を記述した連立方程式を高速に解くために、残差切除法の原理となる反復過程における解の修正量の求め方を流速ベクトル及び圧力を反復的に収束させる演算部分に適用することで、解が収束するまでの時間を飛躍的に減少する非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体を提供することを目的とする。
【００２３】
また本発明は、過去の圧力修正量を線形に結合するパラメータを、求める流速ベクトル及び圧力の厳密解との残差の二乗和の平方根が最小となるよう決定することにより、反復計算での新たな圧力修正量を求める非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体を提供することを目的とする。
【課題を解決するための手段】
【００２４】
上記目的を達成するために第１発明に係る非圧縮性流体の非定常解析装置は、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置において、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶する記憶手段と、前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段とを備えることを特徴とする。
【００２５】
また、第２発明に係る非圧縮性流体の非定常解析装置は、第１発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段とを備えることを特徴とする。
【００２６】
上記目的を達成するために第３発明に係る非圧縮性流体の非定常解析方法は、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求めるコンピュータを用いた非圧縮性流体の非定常解析方法において、前記コンピュータは、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶し、前記解析点での流速ベクトル及び圧力の初期値を受け付け、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出し、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断し、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出し、算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、算出した新たな圧力修正量を算出回数と対応付けて記憶し、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を用いて流速ベクトル近似解を算出し、再度流速ベクトル近似解の発散を算出し、算出した発散の絶対値の最大値が所定値より小さいか否か判断し、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力することを特徴とする。
【００２７】
また、第４発明に係る非圧縮性流体の非定常解析方法は、第３発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出し、算出した残差の二乗和の平方根を最小とするパラメータを算出し、算出したパラメータを用いて、記憶してある圧力修正量を線形に結合することを特徴とする。
【００２８】
上記目的を達成するために第５発明に係るコンピュータプログラムは、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムにおいて、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段とを備えることを特徴とする。
【００２９】
また、第６発明に係るコンピュータプログラムは、第５発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段とを備えることを特徴とする。
【００３０】
上記目的を達成するために第７発明に係るコンピュータプログラムを記録した記録媒体は、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムを記録した記録媒体において、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段とを備えることを特徴とする。
【００３１】
また、第８発明に係るコンピュータプログラムを記録した記録媒体は、第７発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段とを備えることを特徴とする。
【００３２】
第１発明、第３発明、第５発明及び第７発明では、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を求めるべく、単位計算時間毎に、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、圧力の初期値を圧力近似解として、解析点での流速ベクトル近似解を算出し、解析点での流速ベクトル近似解の発散を算出する。解析点にて算出した発散の絶対値の最大値が所定値より大きい場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出し、算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、算出した新たな圧力修正量を算出回数と対応付けて記憶し、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を用いて流速ベクトル近似解を算出し、再度流速ベクトル近似解の発散を算出し、算出した発散の絶対値の最大値が所定値より小さいか否か判断する。流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さい場合、流速ベクトル近似解及び圧力近似解を求める流速ベクトル及び圧力として出力する。これにより、非圧縮性流体の非定常挙動を記述した連立方程式を高速に解くために、残差切除法の原理となる反復過程における解の修正量の求め方を流速ベクトル及び圧力を反復的に収束させる演算部分に適用することができ、反復演算過程において、以前に算出された圧力修正量を用いて、線形にパラメータ結合させることで新たな圧力修正量を求めることにより、流体の連続の式を満足しつつ解が収束するまでの時間を飛躍的に減少することが可能となる。
【００３３】
また、第２発明、第４発明、第６発明及び第８発明では、新たな圧力修正量を求めるために、解析点において、残差として流速ベクトル近似解の発散を単位計算時間で除した値を算出し、算出した残差の二乗和の平方根を最小とするパラメータを算出し、算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する。これにより、新たな圧力修正量を用いた反復演算が発散することなく確実に解を収束させることができ、解を求めるまでの時間を飛躍的に減少することが可能となる。
【発明の効果】
【００３４】
第１発明、第３発明、第５発明又は第７発明によれば、非圧縮性流体の非定常挙動を記述した連立方程式を高速に解くために、残差切除法の原理となる反復過程における解の修正量の求め方を流速ベクトル及び圧力を反復的に収束させる演算部分に適用することができ、反復演算過程において、以前に算出された圧力修正量を用いて、線形にパラメータ結合させることで新たな圧力修正量を求めることにより、流体の連続の式を満足しつつ解が収束するまでの時間を飛躍的に減少することが可能となる。
【００３５】
また、第２発明、第４発明、第６発明又は第８発明によれば、新たな圧力修正量を用いた反復演算が発散することなく確実に解を収束させることができ、解を求めるまでの時間を飛躍的に減少することが可能となる。
【発明を実施するための最良の形態】
【００３６】
以下、本発明の実施の形態について、図面を参照して説明する。図１は、本発明の実施の形態に係る非圧縮性流体の非定常解析装置を具現化するコンピュータの概略構成図である。図１に示すように、非圧縮性流体の非定常解析装置を具現化するコンピュータ１は、少なくとも、ＣＰＵ（中央演算装置）１１、記憶手段１２、ＲＡＭ（メモリ）１３、外部の通信手段と接続する通信手段１４、マウス、キーボード等の入力手段１５、モニタ、プリンタ等の出力手段１６、及び補助記憶手段１７で構成される。
【００３７】
ＣＰＵ１１は、バスを介して上述したようなハードウェア各部と接続されており、上述したハードウェア各部を制御するとともに、記憶手段１２に格納されたプログラムに従って、種々のソフトウェア的機能を実行する。ＲＡＭ（メモリ）１３は、ＳＲＡＭ、フラッシュメモリ等からなり、プログラム処理において一時的に発生したデータを一時記憶する。
【００３８】
出力手段１６は、液晶表示装置、ＣＲＴディスプレイ等の表示装置、又はプリンタ等の印刷装置であり、演算結果を表示し、又は印刷して出力する。
【００３９】
補助記憶手段１７は、非圧縮性流体の非定常解析装置を具現化するコンピュータで使用するプログラムを記録した可搬型記録媒体１８であり、ＤＶＤ、ＣＤ−ＲＯＭ等が該当する。また、非圧縮性流体の非定常解析装置で使用するデータを記録する可搬型記録媒体１８等も含む。
【００４０】
非圧縮性流体の非定常挙動を解析する場合、求めるべき解は、（数８）、（数９）で表されるナビエ・ストークスの運動方程式及び流体の連続の式を同時に満たす流速ベクトル及び圧力である。（数８）、（数９）において、ベクトルｕは流体の流速ベクトルを、ｐは圧力を流体密度で除算した値を、ｔは時間を、νは動粘性係数を、それぞれ示す。
【００４１】
【数８】
2005240889.tif 000010
【００４２】
【数９】
2005240889.tif 000011
【００４３】
図２は、本発明の実施の形態に係る非圧縮性流体の非定常解析装置のＣＰＵ１１での処理のフローチャートである。ＣＰＵ１１は、（数８）、（数９）を解くことにより、流速ベクトルｕ及び圧力ｐを求める。
【００４４】
図２において、ＣＰＵ１１は、流速ベクトルｕ及び圧力ｐの初期値を受け付け（ステップＳ２０１）、流速ベクトルｕの初期近似解を（数１０）に基づいて算出する（ステップＳ２０２）。なお（数１０）において、流速ベクトルｕの上部に付された波線（〜）は流速ベクトルｕの近似解であることを示しており、ｎは所定の時刻における状態であることを示している。
【００４５】
【数１０】
2005240889.tif 000012
【００４６】
ＣＰＵ１１は、（数８）から（数１０）を差し引くことにより、（数１１）を導く。（数１１）では、右辺は所定の時刻から単位計算時間Δｔが経過した後の圧力から、所定の時刻での圧力を引いた圧力修正量δｐの勾配を、左辺は流速ベクトルｕの単位計算時間Δｔ当たりの変化を、それぞれ示している。
【００４７】
【数１１】
2005240889.tif 000013
【００４８】
次に、ＣＰＵ１１は、（数１１）の両辺に対して発散を取り、（数９）を代入することで、（数１２）に示すように、圧力修正量δｐの勾配の発散、すなわち流速ベクトルｕの発散を単位計算時間で除した値を算出する。
【００４９】
【数１２】
2005240889.tif 000014
【００５０】
ｕn+1及びｐn+1を求めるためにＣＰＵ１１は、任意の回数ｋだけ反復演算して圧力修正量Δｐkを求める。すなわち、ＣＰＵ１１は、カウンタｋの初期値を０（ゼロ）とし（ステップＳ２０３）、流速ベクトルｕの初期値を（数３）で算出した流速ベクトルｕの初期近似解とし（ステップＳ２０４）、全ての解析点における流速ベクトルｕkについて発散divｕkの絶対値の最大値が閾値であるεより小さいか否かを判断する（ステップＳ２０５）。ＣＰＵ１１が、流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより大きいと判断した場合（ステップＳ２０５：ＮＯ）、ＣＰＵ１１は、（数１２）の右辺を残差ｒkとし（ステップＳ２０６）、該残差ｒkの二乗和の平方根が最小となる圧力修正量δｐに対するパラメータαl（ｌは１からＬまでの自然数）を算出する（ステップＳ２０７）。なお、残差ｒkの二乗和の平方根は、所定の時刻から単位計算時間Δｔ経過後における残差をｒk+1として、（数１３）に示すように定める。なお、（数１３）において、Ｎは解析点の数を示している。
【００５１】
【数１３】
2005240889.tif 000015
【００５２】
（数１２）の右辺で定義される残差ｒkの二乗和の平方根を最小とするパラメータαl（ｌは１からＬまでの自然数）を定めるには、例えば最小二乗法を適用する。すなわち、（数１４）に示す残差ｒkに対するＬ元連立方程式を解いて、パラメータαlを定めることになる。Ｌの値を十分に大きくした場合、解が発散するおそれは少なくなることから確実にパラメータαlを定めることができるが、演算時間はより長くなる。実際的には、Ｌは４、５程度で十分実用に耐えうる。
【００５３】
【数１４】
2005240889.tif 000016
【００５４】
ＣＰＵ１１は、（数１４）を解くことにより定まった残差ｒkの二乗和の平方根を最小とするパラメータαlを用いて、新たな圧力修正量δｐkを（数１５）により算出する（ステップＳ２０８）。
【００５５】
【数１５】
2005240889.tif 000017
【００５６】
そして、ＣＰＵ１１は、（数１５）で算出した新たな圧力修正量δｐkを用いて、（数１６）により流速ベクトルｕ及び圧力ｐを算出する（ステップＳ２０９）。算出した流速ベクトルｕkを（数３）で算出した流速ベクトルｕの近似解とし（ステップＳ２１０）、カウンタｋを１インクリメントして（ステップＳ２１１）、ステップＳ２０４に戻る。
【００５７】
【数１６】
2005240889.tif 000018
【００５８】
ＣＰＵ１１が、流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより小さいと判断した場合（ステップＳ２０５：ＹＥＳ）、ＣＰＵ１１は、（数１７）に基づいて、所定の時刻から単位計算時間Δｔ経過後の流速ベクトルｕ及び圧力ｐを算出する（ステップＳ２１２）。上述した処理を目的とする時刻ｔmaxまで繰り返すことで、（数８）、（数９）で表されるナビエ・ストークスの運動方程式と流体の連続の式を連立して、所定の時刻における流速ベクトルｕ及び圧力ｐを求めることが可能となる。
【００５９】
【数１７】
2005240889.tif 000019
【００６０】
以上のように、連立方程式を高速に解析できる残差切除法を単純に適用するのではなく、流速ベクトル及び圧力の修正量を決定する演算部分に残差切除法を適用することにより、緩和係数を用いることなく、残差の二乗和の平方根（ノルム）を最小とする反復過程における以前の圧力修正量を線形結合して求める新たな圧力修正量を用いて、効率よく反復処理を実行し、解が収束する時間を短縮する。
【００６１】
以下、具体的な解析例について説明する。図３は本解析に用いた解析格子の例示図であり、図４は本解析に用いた境界条件の例示図である。本解析例では、図３に示すように解析格子の数を６２０００点とし、左側中央に流体を流入させるノズル３１を配している。
【００６２】
境界条件として、ノズル３１の背面４１は、粘性抵抗を有しない壁であると設定し、該背面４１と対向している側４２は、流体が自然流出する面であると設定する。ノズル３１の背面４１は粘性抵抗を有しないことから、流体は背面４１に沿って抵抗を受けずに流れる。
【００６３】
また、流れ方向の側壁面４３、４４では速度０（ゼロ）ｍ／ｓｅｃと設定し、ノズル３１への流入速度を２ｍ／ｓｅｃと設定した。
【００６４】
流速ベクトル及び圧力が共に０（ゼロ）である状態から、所定の時刻に流体を流し始めた場合について、非定常過程の解析を離散化手法として有限体積法を用い、反復解法として、ポイントＳＯＲ法を用いて行った。該解析結果を、従来の解析方法を用いた場合（緩和係数ωを用いる方法）での解析結果と比較する。従来の方法で用いる緩和係数ωは、１から２の範囲の値として設定した。
【００６５】
比較を容易にするため、初期状態から１回目の反復演算終了まで（1ステップ後）、及び１０回目の反復演算終了まで（１０ステップ後）に要した時間を、従来の解析方法を用いた場合、及び本実施の形態に係る解析方法を用いた場合について、それぞれ測定した。図５は、従来の解析方法を用いた場合、及び本実施の形態に係る解析方法を用いた場合における測定結果を示す図である。なお、収束の判定に用いる閾値εとして、十分に小さな値である１０-4を用いた。
【００６６】
図５から明らかなように、解析に要する時間は、本発明に係る解析方法を用いた場合、従来の解析方法を用いた場合と比較して、第１ステップまでに要する時間では、緩和係数ωを最適な値である１．７とした場合であっても約１４倍高速化されていることがわかる。また、１０ステップまでの解析に要する時間で比較した場合、約６．６倍高速化されていることがわかる。追跡するステップ数が大きくなり、流れが定常状態に近づいた場合、解の時間に対する変化が小さくなり、本発明に係る解析方法と従来の解析方法との解析に要する時間の差が小さくなる。しかし、変動が絶えず生じるような非定常挙動を解析する場合、解析に要する時間の差は小さくならず、本実施の形態に係る解析方法の有効性が極めて高いことは明らかである。
【００６７】
また、図５に示すように、従来の解析方法で用いる緩和係数ωを大きく設定した場合、解が収束するまでの時間は短くなっている。しかし、緩和係数ωが２に近づくにつれ収束が不安定となり解が発散している。図５の例では、緩和係数ωを１．８とした場合、解が発散している。このように、緩和係数ωの設定によって、解が収束するか否かが左右されることから、実際には緩和係数ωを最適にして解析することは困難である。
【００６８】
斯かる問題を回避すべく、従来の解析方法では、解が発散しないように緩和係数ωを１近傍に設定することが多く、本実施の形態に係る解析方法との解析速度の差はさらに大きくなる。例えば、緩和係数ωを１とした場合、第１０ステップまでの解析に要する時間の差は、約１０倍となっている。
【００６９】
以上のように本実施の形態によれば、非圧縮性流体の非定常解析に要する時間を従来の解析方法に比べて大幅に短縮することができ、かつ流体の連続の式を満足することを確実に保証することができることから、その意義は極めて高いことは明らかである。
【図面の簡単な説明】
【００７０】
【図１】本発明の実施の形態に係る非圧縮性流体の非定常解析装置を具現化するコンピュータの概略構成図である。
【図２】本発明の実施の形態に係る非圧縮性流体の非定常解析装置のＣＰＵでの処理のフローチャートである。
【図３】本解析に用いた解析格子の例示図である。
【図４】本解析に用いた境界条件の例示図である。
【図５】従来の解析方法を用いた場合、及び本実施の形態に係る解析方法を用いた場合における測定結果を示す図である。
【図６】従来の解析方法に従って、コンピュータを用いて処理した場合の手順を示すフローチャートである。
【符号の説明】
【００７１】
１  コンピュータ
１１  ＣＰＵ
１２  記憶手段
１３  ＲＡＭ
１４  通信手段
１５  入力手段
１６  出力手段
１７  補助記憶手段
１８  可搬型記録媒体
３１  ノズル
【図１】
2005240889.tif 000020
【図２】
2005240889.tif 000021
【図３】
2005240889.tif 000022
【図４】
2005240889.tif 000023
【図５】
2005240889.tif 000024
【図６】
2005240889.tif 000025