JP2005240889A - 非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体 - Google Patents

非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体 Download PDF

Info

Publication number
JP2005240889A
JP2005240889A JP2004050346A JP2004050346A JP2005240889A JP 2005240889 A JP2005240889 A JP 2005240889A JP 2004050346 A JP2004050346 A JP 2004050346A JP 2004050346 A JP2004050346 A JP 2004050346A JP 2005240889 A JP2005240889 A JP 2005240889A
Authority
JP
Japan
Prior art keywords
velocity vector
approximate solution
pressure
flow velocity
calculated
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2004050346A
Other languages
English (en)
Other versions
JP4474943B2 (ja
Inventor
幸司 ▲高▼谷
Koji Takatani
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nippon Steel Corp
Original Assignee
Sumitomo Metal Industries Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sumitomo Metal Industries Ltd filed Critical Sumitomo Metal Industries Ltd
Priority to JP2004050346A priority Critical patent/JP4474943B2/ja
Publication of JP2005240889A publication Critical patent/JP2005240889A/ja
Application granted granted Critical
Publication of JP4474943B2 publication Critical patent/JP4474943B2/ja
Anticipated expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

【課題】 非圧縮性流体の非定常挙動を記述した連立方程式の解が収束するまでの時間を飛躍的に減少する非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体を提供する。
【解決手段】 解析点での流速ベクトル及び圧力の初期値を受け付け、解析点での流速ベクトル近似解を算出し、流速ベクトル近似解の発散の絶対値の最大値が所定値より大きい場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出した圧力修正量と記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいか再度判断し、流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さい場合、流速ベクトル近似解及び圧力近似解を出力する。
【選択図】 図２

Description

本発明は、非圧縮性流体の非定常挙動を有限要素法、差分法、有限体積法等を用いて有限な解析点に離散化して解析する非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体に関する。
従来の非圧縮性流体の非定常挙動を解析する方法について、層流を例に挙げて説明する。なお、例えば乱流であっても、非定常挙動の解析方法として本質的な相違点は生じない。
非圧縮性流体の非定常挙動を解析する場合、解を求めるべき方程式は、（数１）、（数２）で表されるナビエ・ストークスの運動方程式及び流体の連続性条件式である。（数１）、（数２）において、ベクトルｕは流体の流速ベクトルを、ｐは圧力を流体密度で除算した値を、ｔは時間を、νは動粘性係数を、それぞれ示す。
Figure 2005240889
Figure 2005240889
（数１）、（数２）の解を求める手段としては、非特許文献１に開示されているＳＭＡＣ法、非特許文献２に開示されているＦｒａｃｔｉｏｎａｌＳｔｅｐ法等が良く知られている。しかし、いずれの方法においても、解の収束条件に流体の連続条件式が含まれておらず、流体の状況が急変した場合においては、解の精度が保証できないという問題点があった。
そこで、例えば非特許文献３に開示されているＳＯＬＡ法、非特許文献４に開示されているＴＡＫＥＭＩＴＳＵの方法が提案されている。図６は、非特許文献３又は非特許文献４に開示されている従来の解析方法に従って、コンピュータを用いて処理した場合の手順を示すフローチャートである。
図６では、まず流速ベクトルｕ及び圧力ｐの初期値が入力され（ステップＳ６０１）、流速ベクトルｕの初期近似解を（数３）に基づいて算出する（ステップＳ６０２）。なお（数３）において、ベクトルの上部に付された波線（〜）は近似解であることを示しており、ｎは所定の時刻における状態であることを示している。
Figure 2005240889
（数１）から（数３）を差し引くことにより、（数４）を求める。（数４）では、右辺は所定の時刻から単位計算時間Δｔが経過した後の圧力から、所定の時刻での圧力を引いた圧力修正量δｐの勾配、すなわちグラーディエント（gradient）δｐを、左辺は流速ベクトルｕの単位計算時間Δｔ当たりの変化を、それぞれ示している。
Figure 2005240889
次に、（数４）の両辺に対して発散を取り、（数２）を代入して得ることができる（数５）に示すように、圧力修正量δｐの勾配の発散、すなわち流速ベクトルｕの発散を単位時間で除した値を算出する。
Figure 2005240889
ｕn+1及びｐn+1を求めるために、任意の回数ｋだけ反復演算する。すなわち、カウンタｋの初期値を０（ゼロ）とし（ステップＳ６０３）、流速ベクトルｕの初期値を（数３）で算出した流速ベクトルｕの近似解とし（ステップＳ６０４）、全ての解析点における流速ベクトルｕkについて発散divｕkの絶対値の最大値が閾値であるεより小さいか否かを判断する（ステップＳ６０５）。流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより大きいと判断した場合（ステップＳ６０５：ＮＯ）、（数５）より圧力修正量δｐを算出し（ステップＳ６０６）、緩和係数ωを用いて、（数６）のように所定の時刻における流速ベクトルｕk及び圧力ｐkを算出して（ステップＳ６０７）、所定の時刻における流速ベクトルｕkを（数３）で算出した流速ベクトルｕの近似解とし（ステップＳ６０８）、カウンタｋを１インクリメントして（ステップＳ６０９）、ステップＳ６０４へ戻る。
Figure 2005240889
流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより小さいと判断した場合（ステップＳ６０５：ＹＥＳ）、（数７）に示すように、求めた流速ベクトルｕk、圧力ｐkを、単位計算時間Δｔ経過後の流速ベクトルｕn+1、圧力ｐn+1として算出する（ステップＳ６１０）、上述した処理を目的の時刻ｔmaxになるまで繰り返すことで所定の時刻における流体の流速ベクトルｕ及び圧力ｐを求めることができる。
Figure 2005240889
本方法によれば、流体の連続の式を満足しているか否かを、各次官ステップ毎に判定しており（ステップＳ６０５）、全ての時間ステップにおいて、必ず連続の式を満足した解を求めることが可能となる。
エイ．エイ．アムスデン（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における流速ベクトル及び圧力の解を求めるまでに相当の時間を要するという問題点が残されている。
各時間ステップ毎の流速ベクトルと圧力とを求める処理を高速に行うべく、連立方程式の高速解法である残差切除法の原理を利用することも考えられている。最も簡単な方法として、流速ベクトル修正量を求めるステップに残差切除法を適用することが考えられる。
しかし、流速ベクトル修正量は、反復解法を適当な回数だけ反復して求めた近似的な値で十分であり、残差切除法を用いて完全に解を求める場合、たとえ残差切除法を適当な回数で打ち切り近似解を求めるようにした場合であっても、逆に解析時間が増加することは明らかであり好ましくない。
本発明は斯かる事情に鑑みてなされたものであり、非圧縮性流体の非定常挙動を記述した連立方程式を高速に解くために、残差切除法の原理となる反復過程における解の修正量の求め方を流速ベクトル及び圧力を反復的に収束させる演算部分に適用することで、解が収束するまでの時間を飛躍的に減少する非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体を提供することを目的とする。
また本発明は、過去の圧力修正量を線形に結合するパラメータを、求める流速ベクトル及び圧力の厳密解との残差の二乗和の平方根が最小となるよう決定することにより、反復計算での新たな圧力修正量を求める非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体を提供することを目的とする。
上記目的を達成するために第１発明に係る非圧縮性流体の非定常解析装置は、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置において、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶する記憶手段と、前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段とを備えることを特徴とする。
また、第２発明に係る非圧縮性流体の非定常解析装置は、第１発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段とを備えることを特徴とする。
上記目的を達成するために第３発明に係る非圧縮性流体の非定常解析方法は、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求めるコンピュータを用いた非圧縮性流体の非定常解析方法において、前記コンピュータは、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶し、前記解析点での流速ベクトル及び圧力の初期値を受け付け、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出し、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断し、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出し、算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、算出した新たな圧力修正量を算出回数と対応付けて記憶し、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を用いて流速ベクトル近似解を算出し、再度流速ベクトル近似解の発散を算出し、算出した発散の絶対値の最大値が所定値より小さいか否か判断し、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力することを特徴とする。
また、第４発明に係る非圧縮性流体の非定常解析方法は、第３発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出し、算出した残差の二乗和の平方根を最小とするパラメータを算出し、算出したパラメータを用いて、記憶してある圧力修正量を線形に結合することを特徴とする。
上記目的を達成するために第５発明に係るコンピュータプログラムは、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムにおいて、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段とを備えることを特徴とする。
また、第６発明に係るコンピュータプログラムは、第５発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段とを備えることを特徴とする。
上記目的を達成するために第７発明に係るコンピュータプログラムを記録した記録媒体は、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムを記録した記録媒体において、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段とを備えることを特徴とする。
また、第８発明に係るコンピュータプログラムを記録した記録媒体は、第７発明において、前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段とを備えることを特徴とする。
第１発明、第３発明、第５発明及び第７発明では、ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を求めるべく、単位計算時間毎に、求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、圧力の初期値を圧力近似解として、解析点での流速ベクトル近似解を算出し、解析点での流速ベクトル近似解の発散を算出する。解析点にて算出した発散の絶対値の最大値が所定値より大きい場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出し、算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、算出した新たな圧力修正量を算出回数と対応付けて記憶し、算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を用いて流速ベクトル近似解を算出し、再度流速ベクトル近似解の発散を算出し、算出した発散の絶対値の最大値が所定値より小さいか否か判断する。流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さい場合、流速ベクトル近似解及び圧力近似解を求める流速ベクトル及び圧力として出力する。これにより、非圧縮性流体の非定常挙動を記述した連立方程式を高速に解くために、残差切除法の原理となる反復過程における解の修正量の求め方を流速ベクトル及び圧力を反復的に収束させる演算部分に適用することができ、反復演算過程において、以前に算出された圧力修正量を用いて、線形にパラメータ結合させることで新たな圧力修正量を求めることにより、流体の連続の式を満足しつつ解が収束するまでの時間を飛躍的に減少することが可能となる。
また、第２発明、第４発明、第６発明及び第８発明では、新たな圧力修正量を求めるために、解析点において、残差として流速ベクトル近似解の発散を単位計算時間で除した値を算出し、算出した残差の二乗和の平方根を最小とするパラメータを算出し、算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する。これにより、新たな圧力修正量を用いた反復演算が発散することなく確実に解を収束させることができ、解を求めるまでの時間を飛躍的に減少することが可能となる。
第１発明、第３発明、第５発明又は第７発明によれば、非圧縮性流体の非定常挙動を記述した連立方程式を高速に解くために、残差切除法の原理となる反復過程における解の修正量の求め方を流速ベクトル及び圧力を反復的に収束させる演算部分に適用することができ、反復演算過程において、以前に算出された圧力修正量を用いて、線形にパラメータ結合させることで新たな圧力修正量を求めることにより、流体の連続の式を満足しつつ解が収束するまでの時間を飛躍的に減少することが可能となる。
また、第２発明、第４発明、第６発明又は第８発明によれば、新たな圧力修正量を用いた反復演算が発散することなく確実に解を収束させることができ、解を求めるまでの時間を飛躍的に減少することが可能となる。
以下、本発明の実施の形態について、図面を参照して説明する。図１は、本発明の実施の形態に係る非圧縮性流体の非定常解析装置を具現化するコンピュータの概略構成図である。図１に示すように、非圧縮性流体の非定常解析装置を具現化するコンピュータ１は、少なくとも、ＣＰＵ（中央演算装置）１１、記憶手段１２、ＲＡＭ（メモリ）１３、外部の通信手段と接続する通信手段１４、マウス、キーボード等の入力手段１５、モニタ、プリンタ等の出力手段１６、及び補助記憶手段１７で構成される。
ＣＰＵ１１は、バスを介して上述したようなハードウェア各部と接続されており、上述したハードウェア各部を制御するとともに、記憶手段１２に格納されたプログラムに従って、種々のソフトウェア的機能を実行する。ＲＡＭ（メモリ）１３は、ＳＲＡＭ、フラッシュメモリ等からなり、プログラム処理において一時的に発生したデータを一時記憶する。
出力手段１６は、液晶表示装置、ＣＲＴディスプレイ等の表示装置、又はプリンタ等の印刷装置であり、演算結果を表示し、又は印刷して出力する。
補助記憶手段１７は、非圧縮性流体の非定常解析装置を具現化するコンピュータで使用するプログラムを記録した可搬型記録媒体１８であり、ＤＶＤ、ＣＤ−ＲＯＭ等が該当する。また、非圧縮性流体の非定常解析装置で使用するデータを記録する可搬型記録媒体１８等も含む。
非圧縮性流体の非定常挙動を解析する場合、求めるべき解は、（数８）、（数９）で表されるナビエ・ストークスの運動方程式及び流体の連続の式を同時に満たす流速ベクトル及び圧力である。（数８）、（数９）において、ベクトルｕは流体の流速ベクトルを、ｐは圧力を流体密度で除算した値を、ｔは時間を、νは動粘性係数を、それぞれ示す。
Figure 2005240889
Figure 2005240889
図２は、本発明の実施の形態に係る非圧縮性流体の非定常解析装置のＣＰＵ１１での処理のフローチャートである。ＣＰＵ１１は、（数８）、（数９）を解くことにより、流速ベクトルｕ及び圧力ｐを求める。
図２において、ＣＰＵ１１は、流速ベクトルｕ及び圧力ｐの初期値を受け付け（ステップＳ２０１）、流速ベクトルｕの初期近似解を（数１０）に基づいて算出する（ステップＳ２０２）。なお（数１０）において、流速ベクトルｕの上部に付された波線（〜）は流速ベクトルｕの近似解であることを示しており、ｎは所定の時刻における状態であることを示している。
Figure 2005240889
ＣＰＵ１１は、（数８）から（数１０）を差し引くことにより、（数１１）を導く。（数１１）では、右辺は所定の時刻から単位計算時間Δｔが経過した後の圧力から、所定の時刻での圧力を引いた圧力修正量δｐの勾配を、左辺は流速ベクトルｕの単位計算時間Δｔ当たりの変化を、それぞれ示している。
Figure 2005240889
次に、ＣＰＵ１１は、（数１１）の両辺に対して発散を取り、（数９）を代入することで、（数１２）に示すように、圧力修正量δｐの勾配の発散、すなわち流速ベクトルｕの発散を単位計算時間で除した値を算出する。
Figure 2005240889
ｕn+1及びｐn+1を求めるためにＣＰＵ１１は、任意の回数ｋだけ反復演算して圧力修正量Δｐkを求める。すなわち、ＣＰＵ１１は、カウンタｋの初期値を０（ゼロ）とし（ステップＳ２０３）、流速ベクトルｕの初期値を（数３）で算出した流速ベクトルｕの初期近似解とし（ステップＳ２０４）、全ての解析点における流速ベクトルｕkについて発散divｕkの絶対値の最大値が閾値であるεより小さいか否かを判断する（ステップＳ２０５）。ＣＰＵ１１が、流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより大きいと判断した場合（ステップＳ２０５：ＮＯ）、ＣＰＵ１１は、（数１２）の右辺を残差ｒkとし（ステップＳ２０６）、該残差ｒkの二乗和の平方根が最小となる圧力修正量δｐに対するパラメータαl（ｌは１からＬまでの自然数）を算出する（ステップＳ２０７）。なお、残差ｒkの二乗和の平方根は、所定の時刻から単位計算時間Δｔ経過後における残差をｒk+1として、（数１３）に示すように定める。なお、（数１３）において、Ｎは解析点の数を示している。
Figure 2005240889
（数１２）の右辺で定義される残差ｒkの二乗和の平方根を最小とするパラメータαl（ｌは１からＬまでの自然数）を定めるには、例えば最小二乗法を適用する。すなわち、（数１４）に示す残差ｒkに対するＬ元連立方程式を解いて、パラメータαlを定めることになる。Ｌの値を十分に大きくした場合、解が発散するおそれは少なくなることから確実にパラメータαlを定めることができるが、演算時間はより長くなる。実際的には、Ｌは４、５程度で十分実用に耐えうる。
Figure 2005240889
ＣＰＵ１１は、（数１４）を解くことにより定まった残差ｒkの二乗和の平方根を最小とするパラメータαlを用いて、新たな圧力修正量δｐkを（数１５）により算出する（ステップＳ２０８）。
Figure 2005240889
そして、ＣＰＵ１１は、（数１５）で算出した新たな圧力修正量δｐkを用いて、（数１６）により流速ベクトルｕ及び圧力ｐを算出する（ステップＳ２０９）。算出した流速ベクトルｕkを（数３）で算出した流速ベクトルｕの近似解とし（ステップＳ２１０）、カウンタｋを１インクリメントして（ステップＳ２１１）、ステップＳ２０４に戻る。
Figure 2005240889
ＣＰＵ１１が、流速ベクトルｕkの発散divｕkの絶対値の最大値が閾値εより小さいと判断した場合（ステップＳ２０５：ＹＥＳ）、ＣＰＵ１１は、（数１７）に基づいて、所定の時刻から単位計算時間Δｔ経過後の流速ベクトルｕ及び圧力ｐを算出する（ステップＳ２１２）。上述した処理を目的とする時刻ｔmaxまで繰り返すことで、（数８）、（数９）で表されるナビエ・ストークスの運動方程式と流体の連続の式を連立して、所定の時刻における流速ベクトルｕ及び圧力ｐを求めることが可能となる。
Figure 2005240889
以上のように、連立方程式を高速に解析できる残差切除法を単純に適用するのではなく、流速ベクトル及び圧力の修正量を決定する演算部分に残差切除法を適用することにより、緩和係数を用いることなく、残差の二乗和の平方根（ノルム）を最小とする反復過程における以前の圧力修正量を線形結合して求める新たな圧力修正量を用いて、効率よく反復処理を実行し、解が収束する時間を短縮する。
以下、具体的な解析例について説明する。図３は本解析に用いた解析格子の例示図であり、図４は本解析に用いた境界条件の例示図である。本解析例では、図３に示すように解析格子の数を６２０００点とし、左側中央に流体を流入させるノズル３１を配している。
境界条件として、ノズル３１の背面４１は、粘性抵抗を有しない壁であると設定し、該背面４１と対向している側４２は、流体が自然流出する面であると設定する。ノズル３１の背面４１は粘性抵抗を有しないことから、流体は背面４１に沿って抵抗を受けずに流れる。
また、流れ方向の側壁面４３、４４では速度０（ゼロ）ｍ／ｓｅｃと設定し、ノズル３１への流入速度を２ｍ／ｓｅｃと設定した。
流速ベクトル及び圧力が共に０（ゼロ）である状態から、所定の時刻に流体を流し始めた場合について、非定常過程の解析を離散化手法として有限体積法を用い、反復解法として、ポイントＳＯＲ法を用いて行った。該解析結果を、従来の解析方法を用いた場合（緩和係数ωを用いる方法）での解析結果と比較する。従来の方法で用いる緩和係数ωは、１から２の範囲の値として設定した。
比較を容易にするため、初期状態から１回目の反復演算終了まで（1ステップ後）、及び１０回目の反復演算終了まで（１０ステップ後）に要した時間を、従来の解析方法を用いた場合、及び本実施の形態に係る解析方法を用いた場合について、それぞれ測定した。図５は、従来の解析方法を用いた場合、及び本実施の形態に係る解析方法を用いた場合における測定結果を示す図である。なお、収束の判定に用いる閾値εとして、十分に小さな値である１０-4を用いた。
図５から明らかなように、解析に要する時間は、本発明に係る解析方法を用いた場合、従来の解析方法を用いた場合と比較して、第１ステップまでに要する時間では、緩和係数ωを最適な値である１．７とした場合であっても約１４倍高速化されていることがわかる。また、１０ステップまでの解析に要する時間で比較した場合、約６．６倍高速化されていることがわかる。追跡するステップ数が大きくなり、流れが定常状態に近づいた場合、解の時間に対する変化が小さくなり、本発明に係る解析方法と従来の解析方法との解析に要する時間の差が小さくなる。しかし、変動が絶えず生じるような非定常挙動を解析する場合、解析に要する時間の差は小さくならず、本実施の形態に係る解析方法の有効性が極めて高いことは明らかである。
また、図５に示すように、従来の解析方法で用いる緩和係数ωを大きく設定した場合、解が収束するまでの時間は短くなっている。しかし、緩和係数ωが２に近づくにつれ収束が不安定となり解が発散している。図５の例では、緩和係数ωを１．８とした場合、解が発散している。このように、緩和係数ωの設定によって、解が収束するか否かが左右されることから、実際には緩和係数ωを最適にして解析することは困難である。
斯かる問題を回避すべく、従来の解析方法では、解が発散しないように緩和係数ωを１近傍に設定することが多く、本実施の形態に係る解析方法との解析速度の差はさらに大きくなる。例えば、緩和係数ωを１とした場合、第１０ステップまでの解析に要する時間の差は、約１０倍となっている。
以上のように本実施の形態によれば、非圧縮性流体の非定常解析に要する時間を従来の解析方法に比べて大幅に短縮することができ、かつ流体の連続の式を満足することを確実に保証することができることから、その意義は極めて高いことは明らかである。
本発明の実施の形態に係る非圧縮性流体の非定常解析装置を具現化するコンピュータの概略構成図である。 本発明の実施の形態に係る非圧縮性流体の非定常解析装置のＣＰＵでの処理のフローチャートである。 本解析に用いた解析格子の例示図である。 本解析に用いた境界条件の例示図である。 従来の解析方法を用いた場合、及び本実施の形態に係る解析方法を用いた場合における測定結果を示す図である。 従来の解析方法に従って、コンピュータを用いて処理した場合の手順を示すフローチャートである。
符号の説明
１ コンピュータ
１１ ＣＰＵ
１２ 記憶手段
１３ ＲＡＭ
１４ 通信手段
１５ 入力手段
１６ 出力手段
１７ 補助記憶手段
１８ 可搬型記録媒体
３１ ノズル

Claims (8)

  1. ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置において、
    求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶する記憶手段と、
    前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、
    受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、
    前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、
    該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、
    該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、
    算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、
    算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、
    前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段と
    を備えることを特徴とする非圧縮性流体の非定常解析装置。
  2. 前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、
    該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、
    該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段と
    を備えることを特徴とする請求項１に記載の非圧縮性流体の非定常解析装置。
  3. ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求めるコンピュータを用いた非圧縮性流体の非定常解析方法において、
    前記コンピュータは、
    求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶し、
    前記解析点での流速ベクトル及び圧力の初期値を受け付け、
    受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出し、
    前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断し、
    前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出し、
    算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出し、
    算出した新たな圧力修正量を算出回数と対応付けて記憶し、
    算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を用いて流速ベクトル近似解を算出し、再度流速ベクトル近似解の発散を算出し、算出した発散の絶対値の最大値が所定値より小さいか否か判断し、
    前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力することを特徴とする非圧縮性流体の非定常解析方法。
  4. 前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出し、
    算出した残差の二乗和の平方根を最小とするパラメータを算出し、
    算出したパラメータを用いて、記憶してある圧力修正量を線形に結合することを特徴とする請求項３に記載の非圧縮性流体の非定常解析方法。
  5. ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムにおいて、
    求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、
    前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、
    受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、
    前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、
    該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、
    該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、
    算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、
    算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、
    前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段と
    を備えることを特徴とするコンピュータプログラム。
  6. 前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、
    該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、
    該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段と
    を備えることを特徴とする請求項５に記載のコンピュータプログラム。
  7. ナビエ・ストークスの運動方程式及び流体の連続の式を満足する有限個の解析点での流速ベクトル及び圧力を、初期値に基づいて単位計算時間毎に求める非圧縮性流体の非定常解析装置を具現化するコンピュータプログラムを記録した記録媒体において、
    求める圧力と近似圧力との差である圧力修正量を、該圧力修正量を算出する毎に記憶してあり、
    前記解析点での流速ベクトル及び圧力の初期値を受け付ける手段と、
    受け付けた流速ベクトル及び圧力の初期値に基づいて、圧力の初期値を圧力近似解として、前記解析点での流速ベクトル近似解を算出する手段と、
    前記解析点での流速ベクトル近似解の発散を算出し、前記解析点にて算出した発散の絶対値の最大値が所定値より小さいか否か判断する判断手段と、
    該判断手段で、前記流速ベクトル近似解の発散の絶対値の最大値が所定値より大きいと判断した場合、流速ベクトル近似解の発散を単位計算時間で除した値を用いて算出する手段と、
    該手段で算出した圧力修正量と、記憶してある圧力修正量とを線形に結合して新たな圧力修正量を算出する手段と、
    算出した新たな圧力修正量を算出回数と対応付けて記憶する手段と、
    算出した新たな圧力修正量を用いて圧力近似解及び流速ベクトル近似解を算出し、算出した圧力近似解及び流速ベクトル近似解を前記判断手段へ戻す手段と、
    前記判断手段で、前記流速ベクトル近似解の発散の絶対値の絶対値の最大値が所定値より小さいと判断した場合、流速ベクトル近似解及び圧力近似解を出力する手段と
    を備えることを特徴とする記録媒体。
  8. 前記新たな圧力修正量を算出する場合、前記解析点において、流速ベクトル近似解の発散を単位計算時間で除し、残差を算出する手段と、
    該手段で算出した残差の二乗和の平方根を最小とするパラメータを算出する手段と、
    該手段で算出したパラメータを用いて、記憶してある圧力修正量を線形に結合する手段と
    を備えることを特徴とする請求項７に記載の記録媒体。
JP2004050346A 2004-02-25 2004-02-25 非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体 Expired - Fee Related JP4474943B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2004050346A JP4474943B2 (ja) 2004-02-25 2004-02-25 非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2004050346A JP4474943B2 (ja) 2004-02-25 2004-02-25 非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体

Publications (2)

Publication Number Publication Date
JP2005240889A true JP2005240889A (ja) 2005-09-08
JP4474943B2 JP4474943B2 (ja) 2010-06-09

Family

ID=35022838

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2004050346A Expired - Fee Related JP4474943B2 (ja) 2004-02-25 2004-02-25 非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体

Country Status (1)

Country Link
JP (1) JP4474943B2 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021032694A (ja) * 2019-08-23 2021-03-01 正裕 岩永 システム、方法およびプログラム

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021032694A (ja) * 2019-08-23 2021-03-01 正裕 岩永 システム、方法およびプログラム
JP7299613B2 (ja) 2019-08-23 2023-06-28 正裕 岩永 システム、方法およびプログラム

Also Published As

Publication number Publication date
JP4474943B2 (ja) 2010-06-09

Similar Documents

Publication Publication Date Title
Queutey et al. An interface capturing method for free-surface hydrodynamic flows
Abgrall et al. Two-layer shallow water system: a relaxation approach
Pareschi et al. Implicit–explicit Runge–Kutta schemes and applications to hyperbolic systems with relaxation
Caiazzo et al. A numerical investigation of velocity–pressure reduced order models for incompressible flows
Bijl et al. Implicit time integration schemes for the unsteady compressible Navier–Stokes equations: laminar flow
Colmenares et al. Analysis of an augmented mixed‐primal formulation for the stationary B oussinesq problem
Nithiarasu et al. The characteristic‐based split (CBS) scheme—a unified approach to fluid dynamics
Caiazzo Analysis of lattice Boltzmann nodes initialisation in moving boundary problems
Lin et al. A method of lines based on immersed finite elements for parabolic moving interface problems
Barter et al. Shock capturing with higher-order, PDE-based artificial viscosity
Fjordholm et al. Accurate numerical discretizations of non-conservative hyperbolic systems
Li et al. An interfacial lattice Boltzmann flux solver for simulation of multiphase flows at large density ratio
Wong et al. Numerical Stability of Partitioned Approach in Fluid‐Structure Interaction for a Deformable Thin‐Walled Vessel
Murillo et al. Formulation of exactly balanced solvers for blood flow in elastic vessels and their application to collapsed states
Zheng et al. Runge–Kutta–Chebyshev projection method
Murillo et al. Adaptation of flux-based solvers to 2D two-layer shallow flows with variable density including numerical treatment of the loss of hyperbolicity and drying/wetting fronts
Ivančić et al. Energy stable arbitrary Lagrangian Eulerian finite element scheme for simulating flow dynamics of droplets on non–homogeneous surfaces
JP4474943B2 (ja) 非圧縮性流体の非定常解析装置及び方法、コンピュータプログラム、該コンピュータプログラムを記録した記録媒体
Oliver et al. Impact of turbulence model irregularity on high-order discretizations
Zhang et al. A moving mesh finite difference method for non-monotone solutions of non-equilibrium equations in porous media
Eliasson et al. High-order implicit time integration for unsteady turbulent flow simulations
Moro et al. A Hybridized Discontinuous Petrov-Galerkin Method for Compresible Flows
Dzikowski et al. Single component multiphase lattice Boltzmann method for Taylor/Bretherton bubble train flow simulations
Fumagalli et al. Optimal control in ink-jet printing via instantaneous control
Mai-Cao et al. A meshless approach to capturing moving interfaces in passive transport problems

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20060216

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20090623

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20090824

A02 Decision of refusal

Free format text: JAPANESE INTERMEDIATE CODE: A02

Effective date: 20091006

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20091224

A911 Transfer to examiner for re-examination before appeal (zenchi)

Free format text: JAPANESE INTERMEDIATE CODE: A911

Effective date: 20100112

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20100216

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20100301

R150 Certificate of patent or registration of utility model

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130319

Year of fee payment: 3

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130319

Year of fee payment: 3

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313111

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130319

Year of fee payment: 3

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20140319

Year of fee payment: 4

LAPS Cancellation because of no payment of annual fees