(57)【要約】【課題】  流体解析の数値計算方法として、流体の圧縮性の有無や定常・非定常に関係なく、基本方程式を高精度に収束させることができる方法を提供する。【解決手段】  運動量保存則を解き、速度場を近似的に求めるステップ（Ｓ２）と、速度場を用いて各スタッガード点の発散を評価し（Ｓ５）、この発散分布を二次ノルムがより小さくなるようにコントロールして速度修正量を求めるステップ（Ｓ３）と、求めた速度修正量を用いて圧力勾配及び圧力場を修正するステップ（Ｓ４）とを、繰り返し実行し、計算対象領域における速度場および圧力場を求める。【選択図】  図１
2005308735.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータによって、流体解析のための数値計算を行う数値計算方法であって、
計算対象領域に関する初期値と境界条件を受ける第１のステップと、
運動量保存則を解き、前記計算対象領域における速度場を、近似的に求める第２のステップと、
ＭＲ(Mass Relaxation) モジュールによって、前記第２のステップにおいて求めた速度場を用いて各格子点の発散を評価し、この発散分布を二次ノルムがより小さくなるようにコントロールして速度修正量を求める第３のステップと、
前記ＭＲモジュールによって、前記第３のステップにおいて求めた速度修正量を用いて、前記計算対象領域における圧力勾配及び圧力場を修正する第４のステップとを備え、
前記第２〜第４のステップを繰り返し実行し、前記計算対象領域における速度場および圧力場を、求める
ことを特徴とする数値計算方法。
【請求項２】
請求項１において、
前記第３のステップは、
各格子点における発散が一定値となるように、処理するフラット処理と、
前記フラット処理の後、各格子点における発散の二次ノルムが最小になるように、速度修正量を求めるスモール処理とを含む
ことを特徴とする数値計算方法。
【請求項３】
請求項１において、
前記第３のステップは、
各格子点における発散が、零または微小な目標値になるように、速度修正量を求めるδフラット処理を含む
ことを特徴とする数値計算方法。
【請求項４】
請求項１において、
前記第３のステップは、
重み付け加算によって、速度修正量の解の唯一化を行う処理を含む
ことを特徴とする数値計算方法。
【請求項５】
請求項１において、
前記第３のステップは、
スタガード格子点を、発散を評価する格子点として用い、かつ、
各スタガード格子点において求めた速度修正量を基にして、内そうによって、各レギュラー格子点における速度修正量を求める処理と、
内そうによって求めた各レギュラー格子点における速度修正量を、発散の二次ノルムが最小となるように最適化する処理とを含む
ことを特徴とする数値計算方法。
【請求項６】
請求項１において、
前記第４のステップは、
前記第３のステップにおいて求めた速度修正量を用いて、修正速度場を求め、
求めた修正速度場から、修正圧力勾配を求め、
求めた修正圧力勾配を用いて、ノイマンの拘束条件を満たすように、圧力場を修正するものである
ことを特徴とする数値計算方法。
【請求項７】
コンピュータによって、流体解析のための数値計算を行う数値計算装置であって、
計算対象領域に関する初期値と境界条件を受ける第１のステップと、
運動量保存則を解き、前記計算対象領域における速度場を、近似的に求める第２のステップと、
ＭＲ(Mass Relaxation) モジュールによって、前記第２のステップにおいて求めた速度場を用いて各格子点の発散を評価し、この発散分布を二次ノルムがより小さくなるようにコントロールして速度修正量を求める第３のステップと、
前記ＭＲモジュールによって、前記第３のステップにおいて求めた速度修正量を用いて、前記計算対象領域における圧力勾配及び圧力場を修正する第４のステップとを実行するものであり、
前記第２〜第４のステップを繰り返し実行し、前記計算対象領域における速度場および圧力場を、求める
ことを特徴とする数値計算装置。
【請求項８】
コンピュータに、流体解析のための数値計算を行わせる数値計算用プログラムであって、
計算対象領域に関する初期値と境界条件を受ける第１のステップと、
運動量保存則を解き、前記計算対象領域における速度場を、近似的に求める第２のステップと、
ＭＲ(Mass Relaxation) モジュールによって、前記第２のステップにおいて求めた速度場を用いて各格子点の発散を評価し、この発散分布を二次ノルムがより小さくなるようにコントロールして速度修正量を求める第３のステップと、
前記ＭＲモジュールによって、前記第３のステップにおいて求めた速度修正量を用いて、前記計算対象領域における圧力勾配及び圧力場を修正する第４のステップとを備え、
前記第２〜第４のステップを繰り返し実行し、前記計算対象領域における速度場および圧力場を、求める
処理をコンピュータに実行させる数値計算用プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピュータによって、流体解析のための数値計算を行う技術に関するものである。
【背景技術】
【０００２】
流体解析のために用いる基本方程式として、ナビエ・ストークス方程式(Navier-Stokes equation)がある。ナビエ・ストークス方程式は、運動量保存則と連続の式とを含む。そして、これを解く従来の代表的な手法として、非圧縮性の流体については、ＭＡＣ(Marker And Cell)法やＳＩＭＰＬＥ(Semi Implicit Method for Pressure Linked Equation)法があり、また圧縮性のある流体については、ＴＶＤ系の時間進行解法などがある。
【非特許文献１】河村哲也著，「流体解析<２>  数値計算ライブラリ」，朝倉書店，１９９７年５月，pp.49-67
【発明の開示】
【発明が解決しようとする課題】
【０００３】
一般に、非圧縮性の流体の解析においては、運動量保存則と連続の式が、関係式の全てであった。そして、従来技術では、連続の式（発散＝零）を満たすためには、圧力をいかに修正すればよいかが、アルゴリズムの要であった。
【０００４】
本発明は、従来の解法とは全く異なり、流体の圧縮性の有無や、定常・非定常に関係なく実施可能であり、連続の式を含め全ての方程式を高精度に収束させることができる数値計算方法を提供するものである。
【課題を解決するための手段】
【０００５】
本発明は、連続の式すなわち“発散＝零”を用いるのではなく、運動量保存則を解くことによって得られた速度成分を用いて、各点の発散を算定する。そして、計算領域全体の発散の二次ノルムを最小化して、修正速度場を求める。これを用いて修正圧力勾配及び修正圧力場を算定し、運動量保存則に戻る。このプロセスを繰り返して、計算領域全体の発散の二次ノルムを反復毎に小さくし、最終的に、発散＝零  を実現する。本願発明者は、この方法を、質量緩和法（ＭＲ法：Mass Relaxation Method）と呼ぶ。
【０００６】
具体的には本発明は、コンピュータによって、流体解析のための数値計算を行うための方法、装置、およびプログラムとして、計算対象領域に関する初期値と境界条件を受ける第１のステップと、運動量保存則を解き、前記計算対象領域における速度場を、近似的に求める第２のステップと、ＭＲ(Mass Relaxation) モジュールによって、前記第２のステップにおいて求めた速度場を用いて各格子点の発散を評価し、この発散分布を二次ノルムがより小さくなるようにコントロールして速度修正量を求める第３のステップと、前記ＭＲモジュールによって、前記第３のステップにおいて求めた速度修正量を用いて、前記計算対象領域における圧力勾配及び圧力場を修正する第４のステップとを実行し、前記第２〜第４のステップを繰り返し実行し、前記計算対象領域における速度場および圧力場を求めるものである。
【発明の効果】
【０００７】
本発明によると、圧力勾配Ｓm 型運動量保存則を解いて得られた速度成分を用いて、各点の発散を算定し、計算領域全体の発散の二次ノルムを最小化して、修正速度場を求める。これを用いて修正圧力勾配及び修正圧力場を算定し、運動量保存則に戻る。このプロセスを繰り返すことによって、流体解析のための方程式を、高精度に収束させることができる。
【発明を実施するための最良の形態】
【０００８】
以下、本発明の実施形態を図面に基づいて詳細に説明する。
【０００９】
図１は本発明の一実施形態に係る数値計算方法を示すフローチャートである。また、図２は本実施形態に係る数値計算装置を示すブロック図である。本実施形態の数値計算方法は、例えばＣＤ−ＲＯＭ２のような記録媒体に記録されている、当該方法を実現するためのプログラムを、コンピュータに実行させることによって、実施することができる。また、当該方法を実現するためのプログラムを実行するコンピュータを備えた装置によって、容易に実施することができる。
【００１０】
ここでは、説明のために、次のような記号と定義を用いる。
【００１１】
関係式には、物体の境界に沿った境界適合型の一般曲線座標系（ξ1，ξ2，ξ3）を用いることにする。
Notation（全て無次元）
ｘl ：デカルト座標系（ｌ＝１，２，３）
ξM ：一般曲線座標系（Ｍ＝１，２，３）
ｕl ：ｘl軸方向の速度
ＵM ：ξM軸方向の反変速度
ＶM ：ξM軸方向の速度
ｐ  ： 圧力
ρ  ：密度
Ｔ  ：静温度
ｋ  ：乱流エネルギー
μt ：乱流粘性係数
Ｒｅ：レイノルズ数
添字（ｉ，ｊ，ｋ）：ξ1，ξ2，ξ3 軸方向の計算格子を表す。
【００１２】
座標変換のmetricとこれに関係する式は、次のように定義する。
【数１】
2005308735.tif 000003
注  Ｍは総和記号ではないものとする。
【００１３】
図１において、説明を簡単にするために、非圧縮性かつ定常の場合の処理を例にとって、説明を行う。この場合、ステップＳ１１，Ｓ１２，Ｓ１３は省略される。
【００１４】
まずステップＳ１において、計算対象領域に関して適切に設定された初期値および境界条件を受ける。そしてステップＳ２において、圧力勾配Ｓm 型運動量保存則を境界値問題として陰解法によって解き、速度場を近似的に求める。これにより、速度成分の近似解ｕl（ｌ＝１，２，３）が求められる。
【００１５】
次に、ステップＳ５において、ステップＳ２において求めた速度場を用いて、各格子点の発散を評価する。そしてステップＳ３において、発散分布を、その二次ノルムがより小さくなるようにコントロールして、速度修正量を求める。さらにステップＳ４において、ステップＳ３において求めた速度修正量を用いて、圧力勾配及び圧力場を修正する。これらの処理は、数値計算装置１内の演算部１０にあるＭＲモジュール１１によって、実行される。
【００１６】
ステップＳ４において修正された圧力勾配及び圧力場を用いて、ステップＳ２において、再度、運動量保存則を解き、速度場を近似的に求める。このようにステップＳ２〜Ｓ４を、発散の二次ノルムが所定値を下回るまで繰り返し実行する。
【００１７】
各処理の具体的な内容を、数式を交えて説明する。
【００１８】
＜発散の定義と評価＞
本実施形態では、発散を、計算格子であるスタガード格子点(i+1/2,j+1/2,k+1/2)において、次のように定義する。計算領域を、
１≦ｉ≦ＩＭ，１≦ｊ≦ＪＭ，１≦ｋ≦ＫＭ
とし、特記しない限り総和規約を用いることにする。
【数２】
2005308735.tif 000004
【００１９】
一般曲線座標系の反変速度は、メトリックの影響が強く、内そう誤差が大きい。そこで、反変速度の代りにξM軸方向の速度ＶMを導入する。簡単のため、次の定義を用いる。
【数３】
2005308735.tif 000005
但し、Ｍは総和規約ではないものとする。（Ｄ２）を用いれば、（Ｄ１）は次のように書ける。
【数４】
2005308735.tif 000006
ここで、修正前後のξM軸方向の速度場をそれぞれ、ＶM0，ＶMとし、修正量をεMとすると、
【数５】
2005308735.tif 000007
修正前の速度場ＶM0 に対応する発散をＣ0とすると、（Ｄ５）から、
【数６】
2005308735.tif 000008
（Ｄ５）に（Ｄ６）、（Ｄ７）を適用し、Appendix Aを用いれば、
【数７】
2005308735.tif 000009
なお、ＤM（ＶM0 ），ＤM（εM ）の離散化は、appendix Aにおいて、
ＸM＝ＶM0   または  ＸM＝εM
と置いたものである。
【００２０】
【数８】
2005308735.tif 000010
のとき、（Ｄ９）は連続の式を表す。後述するように、既存のアルゴリズムとは異なり、本発明の方法では、反復毎に計算領域の発散の分布を一定にした後で、その二次ノルム‖Ｃ‖が最小になるようにコントロールする。これは、領域内の発散の分布がなだらかで、｜Ｃ｜も小さくなることを意味する。
【００２１】
＜発散のコントロールと速度修正量の決定＞
本実施形態では、後述するフラット処理およびスモール処理によって、発散分布をコントロールし、それに対応する速度修正量を定める方法Ａと、零または微小の目標値δを用いるδフラット処理のみによって、対応する速度修正量を定める方法Ｂとについて、説明する。
【００２２】
（フラット処理：Flat Process）
運動量保存則を解いて得られた速度分布ＶM0 を用いて、（Ｄ７）に従って発散Ｃ0を評価すると、一般には、計算空間において極めて変化の激しい分布となる。このＣ0分布を、後述するスモール処理によって直接ノルムを最小にするのは、必ずしも安定ではない。そこで、発散分布が、計算空間全体において所定の値を持つようにすることを考える。すなわち、次のように発散Ｃの各軸方向に微分が零となるように解けば、発散Ｃは空間全体で一定値とすることができる。
【数９】
2005308735.tif 000011
【００２３】
ξM軸方向の添字をｍとし、その他の添字を省略すれば、ξM軸について（Ｆ１）は離散形で（Ｆ２）のように表される。
【数１０】
2005308735.tif 000012
【００２４】
（Ｄ９）を（Ｆ２）に代入すれば、（Ｆ２）はξ1，ξ2，ξ3 に対し、それぞれ、
【数１１】
2005308735.tif 000013
の三項方程式としてThomas法によって解くことができる。（εM）mのｍ方向の三項方程式の係数はAppendix  Bに与えられている。ε1，ε2，ε3の連立方程式は数回の反復で数値的厳密解に収束する。但し、三つのmatrixを連立した（Ｆ１）は初期値依存性があるので、後述するように唯一化をしなければならない。
【００２５】
得られたスタガード点の解εMを用いると、フラット処理を実行した後の発散は次のように表される。
【数１２】
2005308735.tif 000014
【００２６】
（スモール処理：Small Process）
スモール処理によって、計算領域の発散分布のＬ２ノルム‖Ｃ‖を最小とするように速度修正量ε’M（Ｍ＝１，２，３）を決定する。上述のフラット処理によって、発散ＣFとそのときの修正速度εMが定まっているので、その結果を受けて、現在の発散は次のように表せる。
【数１３】
2005308735.tif 000015
フラット処理を実行したときは、
【数１４】
2005308735.tif 000016
であるが、フラット処理を実行しない場合は、
【数１５】
2005308735.tif 000017
である。 次に、ＣSの二次ノルムを次のように表す。
【数１６】
2005308735.tif 000018
ＣSの二次ノルムを最小にするように修正速度ε’Mを定める。
【数１７】
2005308735.tif 000019
【００２７】
式（Ｓ３）に（Ｄ９）式を代入すると、フラット処理と同様に、ε’1，ε’2，ε’3、に関する三セットの連立方程式の系が得られる。各セットの連立方程式はそれぞれ、ε’1，ε’2，ε’3について三項方程式となり、その係数行列はAppendix Cに示されている。この三セットの連立方程式の系はほんの数回の反復で数値厳密解に到達することができる。但し、フラット処理と同様に、（Ｓ３）式は初期値依存性があるため、解の唯一化をしなければならない。
【００２８】
（δフラット処理）
運動量保存則を解いて得られた速度分布ＶM0を用いて、（D７）に従って発散Ｃ0を評価すると、一般には、計算空間において極めて変化の激しい分布となる。この発散分布が、計算空間全体において目標とする与えられたより小さい一定値δ（零または微小値）となるようにするには、次のように発散Ｃの各軸方向に微分が零となるように解く。その際、不確定境界の境界条件を（Ｆ３）’のようにＣ＝δとすれば、空間全体の発散を、与えられた一定値δとすることが出来る。
【数１８】
2005308735.tif 000020
ξM 軸方向の添字をｍとし、その他の添字を省略すれば、ξM 軸について（Ｆ１）’は離散形で（Ｆ２）’のように表される。
【数１９】
2005308735.tif 000021
【００２９】
不確定境界の境界条件を（Ｆ３）’とする。これは不確定境界方向の修正速度εM の二種条件の形となり、δ＝０のとき、連続の式を満たすことを意味する。
【数２０】
2005308735.tif 000022
（Ｄ９）を（Ｆ２）’に代入し、（Ｆ３）’を適用すれば、（Ｆ２）’はξ1，ξ2，ξ3に対し、それぞれ、
【数１１】
2005308735.tif 000023
の三項方程式としてThomas法によって解く事ができる。（εM）mのｍ方向の三項方程式の係数はAppendix Bに与えられている。ε1，ε2，ε3の連立方程式は数回の反復で数値的厳密解に収束する。但し、三つのMatrixを連立した（Ｆ１）’は初期値依存性があるので、後述するように唯一化をしなければならない。
【００３０】
（最適Staggered grid変換）
Flat & Small処理又はδFlat処理によって得られたStaggered点における修正量εM すなわち（ε1）i j+1/2 k+1/2，（ε2）j k+1/2 i+1/2，（ε3）k i+1/2 j+1/2を適切な方法で内外そうして、格子点における近似値
【数２１】
2005308735.tif 000024
を定める。次に、発散の二次ノルムが最小となるように格子点の近似値
【数２１】
2005308735.tif 000025
を最適化する。近似値の最適化は全点もしくは、｜Ｃ｜の大きい点を選択して実行することもできる。いま、格子点における近似値
【数２１】
2005308735.tif 000026
に重みの（θM ）i j kをつけ、Ｃの二次ノルム‖Ｃ‖2 を最小とするようにθM を定める。Ｈ1 ，Ｈ2 ，Ｈ3 をAppendix F のように定義すると、
【数２２】
2005308735.tif 000027
特定の点（ｉ* ｊ* ｋ* ）に着目すると、全領域のＣの中で
【数２３】
2005308735.tif 000028
を含むのは、点（ｉ* ｊ* ｋ* ）の周りの次の８ブロッグのＣのみである。
【数２４】
2005308735.tif 000029
【数２５】
2005308735.tif 000030
したがって、（G−３）を実行すると、Σの中に微分して残るのは上の８ブロッグのＣのみである。以後、簡単のために、ｉ* ｊ* ｋ* を単にｉｊｋと書くものとすると、（G−３）は（G−４）のように表される。
【数２６】
2005308735.tif 000031
（G−４）の解（θ1）i j k，（θ2）i j k，（θ3）i j kは、（Ｇ−５）式を反復することによって求まる。
【数２７】
2005308735.tif 000032
【００３１】
ここに、
【数２８】
2005308735.tif 000033
【数２９】
2005308735.tif 000034
Ｍ＝１，２，３に対し、それぞれ、ｍ＝ｉ，ｊ，ｋ  ｎ＝ｊ，ｋ，ｉ  ｌ＝ｋ，ｉ，ｊであり、また、複号つきのΣは周辺の４点の和を表す。すなわち、
【数３０】
2005308735.tif 000035
【００３２】
（解の唯一化方法（ＵＡ））
フラット処理およびスモール処理は、ともに、発散を表す微分方程式（Ｄ１）の特性のため、たとえ初期値零からスタートしても、ε1，ε2，ε3またはε’1，ε’2，ε’3の計算順序に応じて、異なる解が存在する。ただし、フラット処理における発散の一定値Ｑ0、およびスモール処理における発散の最小二次ノルムは、いずれも、初期値には依存しない。
【００３３】
さて、いま計算の順序がＳ通りあるとし、そのときの解および重みをそれぞれ、
【数３１】
2005308735.tif 000036
とする。簡単のため、Ｓ＝３とすると、
【数３２】
2005308735.tif 000037
（ＵＡ２）の拘束条件のもとでΩ2が最小となるように重みを定めれば、解の唯一化が達成され、唯一解は、次のように与えられる。
【数３３】
2005308735.tif 000038
【００３４】
（解の唯一化方法（ＵＢ））
Cの二次ノルムは次式で与えられる。
【数３４】
2005308735.tif 000039
【数３５】
2005308735.tif 000040
はStaggered点の
【数３６】
2005308735.tif 000041
をStaggered-gridを変換した回し方の異なるS通りのgrid点の値である。方法UBは（UB−３）のように発散の二次ノルムを最小とする唯一化である。
【数３７】
2005308735.tif 000042
【００３５】
（UB−３）のS元連立方程式から最適
【数３８】
2005308735.tif 000043
が定まり、唯一化された解は（UB−４）で与えられる。
【数３９】
2005308735.tif 000044
【００３６】
（発散コントロールの実際：Ａ）
式（Ａ１）に従って、フラット処理における修正量εMと、スモール処理における修正量ε’Mとを合成し、スタガード格子点における合成修正量を得る。
【数４０】
2005308735.tif 000045
【００３７】
λは加速係数であり、通常は１である。もちろん、フラット処理における修正量εMと、スモール処理における修正量ε’Mとにそれぞれ係数をかけて合成してもかまわない。また、フラット処理を用いない場合は、εM＝０とする。
【００３８】
そして、（Ａ１）を適切な方法で内そうして、格子点（レギュラー格子）における速度修正量
【数４１】
2005308735.tif 000046
を定める。ついで、（Ａ２）に従って、速度場を修正する。
【数４２】
2005308735.tif 000047
【００３９】
なお、ここでは、フラット処理をスモール処理の前処理として実行するものとしたが、スモール処理のみを実行するようにしてもよいし、あるいは、フラット処理のみを実行するようにしてもよい。また、フラット処理とスモール処理を適宜組み合わせて実行してもよい。
【００４０】
（発散コントロールの実際：Ｂ）
δフラット処理によって得られたStaggered点の修正量εM を最適Staggered grid変換によって格子点（レギュラー格子）における速度修正量
【数４１】
2005308735.tif 000048
を定める。ついで（B1）に従って速度場を修正。
【数４３】
2005308735.tif 000049
【００４１】
＜圧力場の修正＞
（Ａ１）または（Ｂ１）の速度ＶM から修正圧力勾配Ｓmを求め、（Ｐ１）から圧力場を修正する。（Ｐ１）はポアソンの方程式で、本発明に係る方法では、境界条件が二種条件となり、典型的なノイマン問題となる。ノイマン問題に対するノイマンの拘束条件（Appendix D参照）を適用して、ポアソンの方程式（Ｐ１）のノイマン問題を解く。
【数４４】
2005308735.tif 000050
ＳmはAppendix Eの圧力勾配型運動量保存則式（１）の離散型のαmlＳmを左辺、それ以外を右辺Ｄlとして得られた次の三元連立一次方程式（Ｐ２）を解いた解（Ｐ３）で、格子点（ｉ，ｊ，ｋ）における値が与えられる。
【数４５】
2005308735.tif 000051
Ｄlはコントロールされた発散分布から得られた修正速度場ＶMから定まるｕ’lを用いて算定する。
【００４２】
なお、ここでは説明を簡単にするために、非圧縮性かつ定常の場合の処理を例にとって、ステップＳ１１，Ｓ１２，Ｓ１３を省略した。ただし、圧縮性流れや、乱流の場合は、ステップＳ１１において、エネルギー式、乱流モデル、状態方程式からそれぞれ、温度Ｔ、乱流粘性係数μt 、乱流エネルギーｋ、および密度ρを算定する。エネルギー式も乱流モデルも、楕円型の境界値問題として陰解法によって解けばよい。また、非定常の場合は、ステップＳ１２，Ｓ１３を実行すればよい。
【００４３】
図３〜図９は本実施形態に係る数値計算方法のコンピュータによる実際の計算結果の例を示すグラフである。図３は計算例の縮流管の形状その他を示す。図４は発散の反復回数による変化を表し、発散の二次ノルムは１０-8まで減少している。図５〜図９は速度及び圧力分布を示している。
【００４４】
＜ＡＰＰＥＮＤＩＸ＞
Appendix A  …  発散の離散化
【数４６】
2005308735.tif 000052
注１：suffixの略記は、積σＸのスタガード点の値はσ及びＸそれぞれのスタガード点の値の積とする。すなわち、
【数４７】
2005308735.tif 000053
を意味する。また、ｉ±１／２，ｊ±１／２，ｋ±１／２はそれぞれ独立に復号同順である。
Appendix B  …  フラット処理およびδフラット処理における三項方程式の係数
【数４８】
2005308735.tif 000054
Appendix C  …  スモール処理の三項方程式の係数
【数４９】
2005308735.tif 000055
Appendix D  …  ノイマンの拘束条件
ａlm＝Ｊｇlm
とおくと、ノイマンの拘束条件は次のようになる。
【数５０】
2005308735.tif 000056
Appendix E  …  圧力勾配Ｓm 型運動量保存側
【数５１】
2005308735.tif 000057
非圧縮の場合は、次のように簡単になる。
【数５２】
2005308735.tif 000058
【００４５】
ｘl 運動量保存則（１）は対流項には風上差分を、その他は中心差分で適切に離散化し、陰解法を適用すれば、定常、非定常に関わりなくｕl に関する楕円型境界値問題となり、圧力勾配Ｓm と他軸方向の速度の近似値及び境界条件を与えれば、速度（ｕl ）ijkを算定することができる。 このmatrixの対角項は、楕円型項のそれと対流項の風上化と非定常なら時間変化項のそれとあいまって優れた優対角性を持つ。 圧力勾配をｐの微分でなく、直接Ｓm で記述したものを、圧力勾配型運動量保存則と呼ぶことにする。
Appendix F・・・・最適化Staggered-grid変換（θM ）i*j*k*を求める方程式の係数ｒ，ＨM ，Ｓp
【数５３】
2005308735.tif 000059
【数５４】
2005308735.tif 000060
次の順に、Ｐ＝１〜７の識別番号を付ける。
(i,j,k),(i+1,j,k),(i-1,j,k),(i,j+1,k),(i,j-1,k),(i,j,k+1),(i,j,k-1)
次の順に、Ｐ＝８〜１９
(i,j+1,k+1),(i,j+1,k-1),(i,j-1,k+1),(i,j-1,k-1)
(i+1,j,k+1),(i+1,j,k-1),(i-1,j,k+1),(i-1,j,k-1)
(i+1,j+1,k),(i+1,j-1,k),(i-1,j+1,k),(i-1,j-1,k)
次の順に、Ｐ＝２０〜２７
(i+1,j+1,k+1),(i+1,j+1,k-1),(i+1,j-1,k+1),(i+1,j-1,k-1)
(i-1,j+1,k+1),(i-1,j+1,k-1),(i-1,j-1,k+1),(i-1,j-1,k-1)
ＳP ：Ｐで表されるsuffixの点に接したstaggered点のｒ2 の和をＳP とする。すなわち、
【数５５】
2005308735.tif 000061
【数５６】
2005308735.tif 000062
【数５７】
2005308735.tif 000063
【産業上の利用可能性】
【００４６】
本発明では、流体解析のための数値計算において、連続の式を含めた基本方程式を高精度に収束させることができるので、例えば、流体解析の数値計算の精度向上や速度向上に、有効である。
【図面の簡単な説明】
【００４７】
【図１】本発明の一実施形態に係る数値計算方法を示すフローチャートである。
【図２】本発明の一実施形態に係る数値計算方法を示すブロック図である。
【図３】数値計算の一例として用いた縮流管の形状等を示す図である。
【図４】数値計算の結果の一例であり、発散の二次ノルムのグラフ表示である。
【図５】数値計算の結果の一例であり、主流速度分布のグラフ表示である。
【図６】数値計算の結果の一例であり、圧力分布のグラフ表示である。
【図７】数値計算の結果の一例であり、圧力分布のグラフ表示である。
【図８】数値計算の結果の一例であり、断面圧力分布の等高線表示である。
【図９】数値計算の結果の一例であり、子午面速度分布のベクトル表示である。
【符号の説明】
【００４８】
１  数値計算装置
１０  演算部
１１  ＭＲモジュール
【図１】
2005308735.tif 000064
【図２】
2005308735.tif 000065
【図３】
2005308735.tif 000066
【図４】
2005308735.tif 000067
【図５】
2005308735.tif 000068
【図６】
2005308735.tif 000069
【図７】
2005308735.tif 000070
【図８】
2005308735.tif 000071
【図９】
2005308735.tif 000072