JP5268496B2 - 流動解析方法、流動解析装置、及び流動解析プログラム - Google Patents

Description

本発明は、解析領域内でスケールが大きく変化する非圧縮流体の挙動解析に適用して好適な、流動解析方法、流動解析装置、及び流動解析プログラムに関する。

従来より、流動解析方法の一つとして、均一な大きさの複数の粒子により非圧縮流体を表現し、複数の粒子の挙動をコンピュータを用いて計算することにより非圧縮流体の挙動を解析する方法、いわゆるMPS（Moving Particle Semi-implicit）法が知られている。MPS法は、空気や水等の自由表面を有する物体、特に大きく変形したり、分裂，合体したりする物体の挙動解析に適しており、原子力分野を始めとしてプラント，航空，船舶，環境等の様々な技術分野において広く利用されている。
特開平７−３３４４８４号公報 特開２０００−３３９２９２号公報 特開２００２−１３７２７２号公報

従来のMPS法では、解析領域内で粒子の大きさ、換言すれば空間解像度を変化させることができない。このため、例えば非常に小さな溝がある流路中における非圧縮流体の挙動を解析する場合、溝のスケールに合わせて粒子の大きさを設定し、設定した均一の大きさの粒子を解析領域内に配置しなければならない。このような背景から、従来のMPS法により解析領域内でスケールが大きく変化する非圧縮流体の挙動を解析する場合には、解析領域内に非常に多くの粒子（計算点）が配置されるために、コンピュータの処理能力の制約上、現実的な時間で流動解析を行うことが困難になる。

本発明は、このような課題を解決するためになされたものであり、その目的は、解析領域内でスケールが大きく変化する非圧縮流体の挙動解析に要する時間を短縮可能な、流動解析方法、流動解析装置、及び流動解析プログラムを提供することにある。

本発明に係る流動解析方法は、非圧縮流体の挙動を解析するための流動解析方法であって、粒子法により定義される粒子を利用して前記非圧縮流体を粒子径が異なる複数の粒子により表現し、かつ、速度の発散が０であることを非圧縮条件とし、Ｖｉ及びＶｊをそれぞれ粒子ｉ及び粒子ｊの体積、ｒｉ及びｒｊをそれぞれ粒子ｉ及び粒子ｊの位置座標、Ｒｉ及びＲｊをそれぞれ粒子ｉ及び粒子ｊの影響半径とした時、以下の数式１で表される粒子ｉと粒子ｊ間の重みｗｉｊを用いて非圧縮流体の挙動を解析することを特徴とする。
本発明に係る流動解析装置は、非圧縮流体の挙動を解析するための流動解析装置であって、粒子法により定義される粒子を利用して前記非圧縮流体を粒子径が異なる複数の粒子により表現する初期データ設定手段を備え、前記非圧縮流体を粒子径が異なる複数の粒子により表現し、かつ、速度の発散が０であることを非圧縮条件とし、Ｖｉ及びＶｊをそれぞれ粒子ｉ及び粒子ｊの体積、ｒｉ及びｒｊをそれぞれ粒子ｉ及び粒子ｊの位置座標、Ｒｉ及びＲｊをそれぞれ粒子ｉ及び粒子ｊの影響半径とした時、以下の数式１で表される粒子ｉと粒子ｊ間の重みｗｉｊを用いて非圧縮流体の挙動を解析することを特徴とする。
本発明に係る流動解析プログラムは、非圧縮流体の挙動を解析するための流動解析プログラムであって、粒子法により定義される粒子を利用して前記非圧縮流体を粒子径が異なる複数の粒子により表現する処理をコンピュータに実行させるに際し、前記非圧縮流体を粒子径が異なる複数の粒子により表現し、かつ、速度の発散が０であることを非圧縮条件とし、Ｖｉ及びＶｊをそれぞれ粒子ｉ及び粒子ｊの体積、ｒｉ及びｒｊをそれぞれ粒子ｉ及び粒子ｊの位置座標、Ｒｉ及びＲｊをそれぞれ粒子ｉ及び粒子ｊの影響半径とした時、以下の数式１で表される粒子ｉと粒子ｊ間の重みｗｉｊを用いることを特徴とする。

本発明に係る流動解析方法、流動解析装置、及び流動解析プログラムによれば、解像度が必要な領域内における粒子の粒子径を小さくする一方で解像度が必要でない領域内における粒子の粒子径を大きくすることにより解析領域内で解像度を変化させることができるので、解析領域内でスケールが大きく変化する非圧縮流体の挙動解析に要する時間を大幅に短縮することができる。

以下、図面を参照して、本発明の実施形態となる流動解析装置の構成及びその動作について説明する。

〔流動解析装置の構成〕
本発明の実施形態となる流動解析装置１は、パーソナルコンピュータ，ワークステーション，汎用コンピュータ等の公知の情報処理装置により構成され、図１に示すように、ＲＡＭ（Random Access Memory）２，ＲＯＭ（Read Only Memory）３，初期設定データベース（ＤＢ）４，及びＣＰＵ（Central Processing Unit）５を主な構成要素として備える。ＲＡＭ２は、ＣＰＵ５が実行する流動解析装置１の起動プログラムや流動解析プログラム等のコンピュータプログラム及び各種データを一時的に格納するワーキングエリアを提供する。ＲＯＭ３は、コンピュータプログラム及びこのコンピュータプログラムを実行するために必要な各種データを記憶する。なおＲＯＭ３は、磁気的又は光学的な記憶媒体若しくは半導体メモリ等のＣＰＵ５により読み取り可能な記憶媒体を含む構成を有し、ＲＯＭ３に記憶されているコンピュータプログラムや各種データの一部又は全部は電気通信回線を介してダウンロードされるように構成してもよい。

初期設定ＤＢ４は、MPS法により流動解析を実行するために必要な各種パラメータを格納する。具体的には、このパラメータには、MPS法における粒子ｉの粒子径ｒ_iの初期値ｒ_iini，位置（重心の座標），速度，温度，圧力，流体の物性（例えば粘度，密度，熱伝導率，比熱、潜熱等），時間刻み幅Δｔ，計算終了時間Ｔ等の初期条件に関する情報が含まれる。ＣＰＵ５は、ＲＯＭ３に記憶されているコンピュータプログラムをＲＡＭ２内にロードし、読み出されたコンピュータプログラムを実行することにより、流動解析装置１全体の処理動作を制御する。具体的には、ＣＰＵ５は、ＲＯＭ３に記憶された流動解析プログラム及び流動解析プログラムを実行するために必要なデータをＲＡＭ２へロードし、流動解析プログラムに従ってMPS法を利用した流動解析を実行する。本実施形態では、ＣＰＵ５は、ＲＯＭ３に記憶されている流動解析プログラムを実行することにより、外力項計算部１１，粘性項計算部１２，粒子移動追跡部１３，圧力項計算部１４，粒子データ更新部１５，及び終了判定部１６として機能する。各部の機能の詳細については後述する。

〔流動解析方法〕
このような構成を有する流動解析装置１は、以下に示す流動解析処理を実行することによりスケールが大きく変化する解析領域内における非圧縮流体の挙動を現実的な時間で解析する。以下、図２に示すフローチャートを参照して、この流動解析処理を実行する際の流動解析装置１の動作について詳しく説明する。

図２に示すフローチャートは、流動解析装置１に対しオペレータが流動解析処理の開始を指示したタイミングで開始となり、流動解析処理はステップＳ１の処理に進む。

ステップＳ１の処理では、ＣＰＵ５が、初期設定ＤＢ４からＲＡＭ２内に解析対象となる非圧縮流体の物性，時間刻み幅Δｔ，計算終了時間Ｔ等の計算条件パラメータ、及び粒子の粒子径，初期座標位置，速度，温度等の初期条件パラメータを読み出し、読み出されたパラメータに基づいて非圧縮流体を粒子径が異なる複数の粒子により表現する。具体的には、ＣＰＵ５は、解像度が必要な領域内には小さな粒子径の粒子を配置し、解像度が必要でない領域内には大きな粒子径の粒子を配置する。これにより、ステップＳ１の処理は完了し、流動解析処理はステップＳ２の処理に進む。

ステップＳ２の処理では、ＣＰＵ５が、現在の座標位置において粒子が適切な粒子径となるように粒子を合体又は分裂させることにより、非圧縮流体の解像度を空間的に変化させる。具体的には、粒子ｉの現在の粒子径Ｌｉに対してその座標位置における理想粒子径をＬ＊とすると、粒子径Ｌｉが理想粒子径Ｌ＊よりも例えば２倍程度大きい場合、ＣＰＵ５は、粒子ｉを２つに分裂させることにより粒子ｉの粒子径Ｌｉを小さくする。一方、粒子径Ｌｉが理想粒子径Ｌ＊よりも例えば２倍程度小さい場合には、ＣＰＵ５は粒子ｉとその最近傍にある粒子ｊを合体させることにより粒子ｉの粒子径Ｌｉを大きくする。なお、この処理により新たに生成される粒子径は粒子の体積を保存させるように決定する。これにより、ステップＳ２の処理は完了し、流動解析処理はステップＳ３の処理に進む。

ステップＳ３の処理では、外力項計算部１１が、重力等の各粒子に直接的に作用する外力項を計算する。粒子法における外力項の計算方法は本願発明の出願時点で既に公知であるので詳細な説明は省略する。これにより、ステップＳ３の処理は完了し、流動解析処理はステップＳ４の処理に進む。

ステップＳ４の処理では、粘性項計算部１２が、粘性項を計算し、粒子ｉが有する速度分布を近傍の粒子ｊに分配することにより、粘性による粒子の速度の拡散を計算する。具体的には、粘性項は以下の数式２により表される。数式２中、パラメータｕ，ｔ，υはそれぞれ速度ベクトル，計算時刻，動粘性係数を示す。この数式２を粒子法で離散化すると、粘性項は陽的には数式３，陰的には数式４のように表される。数式３，４中、パラメータｄ，λ，ｎ^０，ｗ_ｉｊ，ｋはそれぞれ次元数，拡散を解析解に一致させるための係数，基準粒子密度，粒子ｉ，ｊ間の重み，及びタイムステップを示す。

そして粒子ｉの速度変化Δｕは、以下の数式５により表されるので、これを数式４に代入することにより、以下の数式６に示す連立一次方程式が導かれる。本実施形態では、この数式６を解くことにより粘性項を陰的に解く。このように粘性項を陰的に解くことにより、タイムステップｋを大きくすることができるので、発散することなく計算を安定的に行うことができる。

なお本実施形態では、粒子ｉ，ｊ間の重みｗ_ｉｊは、従来のMPS法における重みとは異なり、以下の数式７のように表される。なお数式７中、Ｖ_ｉ及びＶ_ｊはそれぞれ粒子ｉ及び粒子ｊの体積、ｒ_ｉ及びｒ_ｊはそれぞれ粒子ｉ及び粒子ｊの座標、Ｒ_ｉ及びＲ_ｊはそれぞれ粒子ｉ及び粒子ｊの影響半径を示す。これにより、ステップＳ４の処理は完了し、流動解析処理はステップＳ５の処理に進む。

ステップＳ５の処理では、粒子移動追跡部１３が、ステップＳ３及びステップＳ４の計算結果に基づいて粒子の速度を仮更新する。これにより、ステップＳ５の処理は完了し、流動解析処理はステップＳ６の処理に進む。

ステップＳ６の処理では、粒子移動追跡部１３が、ステップＳ５の処理により仮更新された速度による時間刻み幅Δｔの間での粒子の移動を追跡し、追跡結果に基づいて粒子の座標位置を仮更新する。これにより、ステップＳ６の処理は完了し、流動解析処理はステップＳ７の処理に進む。

ステップＳ７の処理では、圧力項計算部１４が、ステップＳ６の処理により算出された座標位置における各粒子の粒子数密度を算出し、算出された粒子数密度が所定値であるか否かを判別する。判別の結果、粒子数密度が所定値でない場合、圧力項計算部１４は、粒子数密度が所定値になるように粒子の圧力を修正することにより粒子の座標位置及び流速を修正する。そして圧力計算部１４は、修正された圧力によって生じる圧力項を算出する。具体的には、圧力項計算部１４は、以下の数式８で表される圧力のポアソン方程式を用いて圧力項を算出する。一般的な圧力のポアソン方程式は、粒子数密度が一定であるとの非圧縮条件に基づいて基準粒子数密度を用いて表現されている。しかしながら本実施形態では、圧力のポアソン方程式は、数式８に示すように、非圧縮状態であっても粒子数密度が一定でない、換言すれば、速度の発散が０であるという非圧縮条件に基づいて基準粒子数密度を用いずに表現されている。これにより、ステップＳ７の処理は完了し、流動解析処理はステップＳ８の処理に進む。

ステップＳ８の処理では、粒子データ更新部１５が、ステップＳ７の処理の算出結果に基づいてＲＡＭ２内に格納されている粒子の速度を更新する。これにより、ステップＳ８の処理は完了し、流動解析処理はステップＳ９の処理に進む。

ステップＳ９の処理では、粒子データ更新部１５が、ステップＳ７の処理の算出結果に基づいてＲＡＭ２内に格納されている粒子の座標位置を更新する。これにより、ステップＳ９の処理は完了し、流動解析処理はステップＳ１０の処理に進む。

ステップＳ１０の処理では、終了判定部１６が、解析時間ｔを時間刻み幅Δｔだけインクリメントする。これにより、ステップＳ１０の処理は完了し、流動解析処理はステップＳ１１の処理に進む。

ステップＳ１１の処理では、終了判定部１６が、解析時間ｔが計算終了時刻Ｔであるか否かを判別する。判別の結果、解析時間ｔが計算終了時刻Ｔである場合、終了判定部１０は、ＲＡＭ２内に格納されている粒子データをＲＯＭ３の所定ファイルにコピーし、一連の流動解析処理を終了する。一方、解析時間ｔが計算終了時刻Ｔでない場合には、終了判定部１０は流動解析処理をステップＳ２の処理に戻す。

〔まとめ〕
本願発明の発明者は、従来のMPS法を利用した流動解析手法において粒子径が異なる複数の粒子により非圧縮流体を表現した場合、図３に示すように圧力値が空間的に振動することを確認した。そして鋭意研究を重ねてきた結果、数式７に示す重みを用いた場合には、図４に示すように滑らかな圧力分布が得られることを知見した。また本願発明の発明者は、粒子間で粒子径が１．５倍以上異なる場合には、図５（ａ）〜（ｄ）に示すように計算が不安定になることを確認した。そして鋭意研究を重ねてきた結果、速度の発散が０であるとの非圧縮条件を与えて流動解析を行うことにより、図６（ａ）〜（ｃ）に示すように計算の安定性が劇的に向上することを知見した。以上のことから、本願発明の発明者は、数式７に示す重みを用いて、速度の発散が０であるとの非圧縮条件を与えて流動解析を行うことにより、解析領域内で解像度を変化させることができ、結果として、スケールが大きく変化する非圧縮流体の挙動解析に要する時間を大幅に短縮することができた。

〔実験例〕
最後に、ダム崩壊問題の解析により本手法の優位性を示す。図７は、４種類の大きさの粒子が混在しているモデルである。粒子の種類を２種類から４種類まで変化させて解析を行ったところ、必要とした粒子数及び計算時間は図８に示すようになった。図８から明らかなように、粒子径を大きく変化させる程、計算効率が高くなり、計算時間は１／３以下となった。また容器右下の点における圧力値の変化を測定した所、図９に示すように粒子径が均一である場合（均一粒子）と粒子径が均一でない場合（非均一粒子）とではほぼ同じ結果になることが確認された。以上のことから、本実施形態の流動解析処理によれば、必要な部分の精度を保ったまま計算コストを下げることが可能であるといえる。

以上、本発明者によってなされた発明を適用した実施の形態について説明したが、この実施形態による本発明の開示の一部をなす記述及び図面により本発明は限定されることはない。例えば、本実施形態における流動解析装置１は、モデム等の通信用部材を備えてインターネットを含む電気通信回線と接続可能なように構成してもよい。この場合には、本実施形態における流動解析プログラムは、電気通信回線からダウンロード等によって、ＲＡＭ２又はＲＯＭ３内に格納させてもよい。なお、この場合における通信ネットワークからダウンロードするためのダウンロードプログラムは、予めＲＯＭ３内に格納されているか、別の記憶媒体からＲＯＭ３内にインストールされるものとする。このように、本実施形態に基づいて当業者等によりなされる他の実施の形態、実施例及び運用技術等は全て本発明の範疇に含まれる。

本発明は、解析領域内でスケールが大きく変化する非圧縮流体の挙動解析に適用することができる。

本発明の実施形態となる流動解析装置の構成を示すブロック図である。 本発明の実施形態となる流動解析処理の流れを示すフローチャート図である。 従来の流動解析処理を用いて粒子径が異なる複数の粒子により表現された非圧縮流体の圧力分布を解析した例を示す図である。 本発明の実施形態となる流動解析処理を用いて粒子径が異なる複数の粒子により表現された非圧縮流体の圧力分布を解析した例を示す図である。 速度の発散が０であるとの非圧縮条件を与えない場合の流動解析例を示す図である。 速度の発散が０であるとの非圧縮条件を与えた場合の流動解析例を示す図である。 ４種類の大きさの粒子が混在している流動解析モデルの一例を示す図である。 図７に示す流動解析モデルについて粒子径を均一にして流動解析を行った場合と不均一にして流動解析を行った場合における粒子数及び計算時間を示す図である。 図７に示す流動解析モデルについて粒子径を均一にして流動解析を行った場合と不均一にして流動解析を行った場合における圧力値の時間変化を示す図である。

符号の説明

１…流動解析装置，２…ＲＡＭ，３…ＲＯＭ，４…初期設定ＤＢ，５…ＣＰＵ，１１…外力項計算部，１２…粘性項計算部，１３…粒子移動追跡部，１４…圧力項計算部，１５…粒子データ更新部，１６…終了判定部

Claims (9)

1. 非圧縮流体の挙動を解析するための流動解析方法であって、
   粒子法により定義される粒子を利用して前記非圧縮流体を粒子径が異なる複数の粒子により表現し、かつ、速度の発散が０であることを非圧縮条件とし、
   Ｖｉ及びＶｊをそれぞれ粒子ｉ及び粒子ｊの体積、ｒｉ及びｒｊをそれぞれ粒子ｉ及び粒子ｊの位置座標、Ｒｉ及びＲｊをそれぞれ粒子ｉ及び粒子ｊの影響半径とした時、以下の数式１で表される粒子ｉと粒子ｊ間の重みｗｉｊを用いて非圧縮流体の挙動を解析することを特徴とする流動解析方法。
   
2. 請求項１に記載の流動解析方法において、前記複数の粒子を合体及び分裂させることにより非圧縮流体の解像度を空間的に変化させることを特徴とする流動解析方法。
3. 請求項１または２に記載の流動解析方法において、前記粒子としてMPS法により定義される粒子を利用することを特徴とする流動解析方法。
4. 非圧縮流体の挙動を解析するための流動解析装置であって、
   粒子法により定義される粒子を利用して前記非圧縮流体を粒子径が異なる複数の粒子により表現する初期データ設定手段を備え、
   前記非圧縮流体を粒子径が異なる複数の粒子により表現し、かつ、速度の発散が０であることを非圧縮条件とし、
   Ｖｉ及びＶｊをそれぞれ粒子ｉ及び粒子ｊの体積、ｒｉ及びｒｊをそれぞれ粒子ｉ及び粒子ｊの位置座標、Ｒｉ及びＲｊをそれぞれ粒子ｉ及び粒子ｊの影響半径とした時、以下の数式１で表される粒子ｉと粒子ｊ間の重みｗｉｊを用いて非圧縮流体の挙動を解析することを特徴とする流動解析装置。
   
5. 請求項４に記載の流動解析装置において、前記複数の粒子を合体及び分裂させることにより非圧縮流体の解像度を空間的に変化させることを特徴とする流動解析装置。
6. 請求項４または５に記載の流動解析装置において、前記粒子としてMPS法により定義される粒子を利用することを特徴とする流動解析装置。
7. 非圧縮流体の挙動を解析するための流動解析プログラムであって、
   粒子法により定義される粒子を利用して前記非圧縮流体を粒子径が異なる複数の粒子により表現する処理をコンピュータに実行させるに際し、
   前記非圧縮流体を粒子径が異なる複数の粒子により表現し、かつ、速度の発散が０であることを非圧縮条件とし、
   Ｖｉ及びＶｊをそれぞれ粒子ｉ及び粒子ｊの体積、ｒｉ及びｒｊをそれぞれ粒子ｉ及び粒子ｊの位置座標、Ｒｉ及びＲｊをそれぞれ粒子ｉ及び粒子ｊの影響半径とした時、以下の数式１で表される粒子ｉと粒子ｊ間の重みｗｉｊを用いることを特徴とする流動解析プログラム。
   
8. 請求項７に記載の流動解析プログラムにおいて、前記複数の粒子を合体及び分裂させることにより非圧縮流体の解像度を空間的に変化させることを特徴とする流動解析プログラム。
9. 請求項７または８に記載の流動解析プログラムにおいて、前記粒子としてMPS法により定義される粒子を利用することを特徴とする流動解析プログラム。