ホーム > 特許ランキング > 住友重機械工業株式会社
知財求人 - 知財ポータルサイト「IP Force」
特開2024-80134シミュレーション方法、シミュレーション装置、及びプログラム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024080134
(43)【公開日】2024-06-13
(54)【発明の名称】シミュレーション方法、シミュレーション装置、及びプログラム
(51)【国際特許分類】
G06Q 99/00 20060101AFI20240606BHJP
【ＦＩ】
G06Q99/00
【審査請求】未請求
【請求項の数】10
【出願形態】ＯＬ
(21)【出願番号】P 2022193059
(22)【出願日】2022-12-01
(71)【出願人】
【識別番号】000002107
【氏名又は名称】住友重機械工業株式会社
(74)【代理人】
【識別番号】100105887
【弁理士】
【氏名又は名称】来山 幹雄
(72)【発明者】
【氏名】廣田 真人
【テーマコード（参考）】
5L049
5L050
【Ｆターム（参考）】
5L049DD02
5L050DD02
(57)【要約】
【課題】混相流体の挙動を安定に解析することが可能なシミュレーション方法を提供する。
【解決手段】混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させる。複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、対象粒子に対する参照粒子の相対速度を小さくして計算する。
【選択図】図３
【解決手段】混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させる。複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、対象粒子に対する参照粒子の相対速度を小さくして計算する。
【選択図】図３
【特許請求の範囲】
【請求項1】
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーション方法であって、
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するシミュレーション方法。
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するシミュレーション方法。
【請求項2】
前記運動量保存式は、解析における数値的な不安定性を低減させる人工粘性項を含んでおり、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、前記対象粒子に付与された密度と、前記複数の参照粒子のそれぞれに付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える請求項１に記載のシミュレーション方法。
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、前記対象粒子に付与された密度と、前記複数の参照粒子のそれぞれに付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える請求項１に記載のシミュレーション方法。
【請求項3】
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーション方法であって、
前記運動量保存式は、解析における数値的な不安定性を低減させる人工粘性項を含んでおり、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるシミュレーション方法。
前記運動量保存式は、解析における数値的な不安定性を低減させる人工粘性項を含んでおり、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるシミュレーション方法。
【請求項4】
前記混相流体に粉粒体が混在しており、前記複数の計算粒子は、前記粉粒体を表す粒子を含んでおり、
前記混相流体を表す前記複数の計算粒子と、前記粉粒体を表す前記複数の計算粒子との挙動を連成解析する請求項１乃至３のいずれか１項に記載のシミュレーション方法。
前記混相流体を表す前記複数の計算粒子と、前記粉粒体を表す前記複数の計算粒子との挙動を連成解析する請求項１乃至３のいずれか１項に記載のシミュレーション方法。
【請求項5】
粒子法を用いて混相流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算する機能と
を有するシミュレーション装置。
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算する機能と
を有するシミュレーション装置。
【請求項6】
前記運動量保存式は、解析における数値的な不安定性を低減させる人工粘性項を含んでおり、
前記処理部は、前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、前記対象粒子の密度と、前記複数の参照粒子の密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える機能を、さらに有する請求項５に記載のシミュレーション装置。
前記処理部は、前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、前記対象粒子の密度と、前記複数の参照粒子の密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える機能を、さらに有する請求項５に記載のシミュレーション装置。
【請求項7】
粒子法を用いて流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式、及び解析における数値的な不安定性を低減させる人工粘性項を含む運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす複数の計算粒子である参照粒子のそれぞれに付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える機能と
を有するシミュレーション装置。
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式、及び解析における数値的な不安定性を低減させる人工粘性項を含む運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす複数の計算粒子である参照粒子のそれぞれに付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える機能と
を有するシミュレーション装置。
【請求項8】
前記混相流体に粉粒体が混在しており、前記複数の計算粒子は、前記粉粒体を表す粒子を含んでおり、
前記処理部は、前記混相流体を表す前記複数の計算粒子と、前記粉粒体を表す前記複数の計算粒子との挙動を連成解析する機能を、さらに有する請求項５乃至７のいずれか１項に記載のシミュレーション装置。
前記処理部は、前記混相流体を表す前記複数の計算粒子と、前記粉粒体を表す前記複数の計算粒子との挙動を連成解析する機能を、さらに有する請求項５乃至７のいずれか１項に記載のシミュレーション装置。
【請求項9】
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーションの手順をコンピュータに実行させるプログラムであって、
前記複数の計算粒子のそれぞれについて質量保存式を解く手順において、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するプログラム。
前記複数の計算粒子のそれぞれについて質量保存式を解く手順において、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するプログラム。
【請求項10】
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式、及び解析における数値的な不安定性を低減させる人工粘性項を含む運動量保存式を含む支配方程式を解いて時間発展させるシミュレーションの手順をコンピュータに実行させるプログラムであって、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く手順において、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるプログラム。
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く手順において、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、シミュレーション方法、シミュレーション装置、及びプログラムに関する。
【背景技術】
【0002】
流体の運動を粒子系の運動として近似して、流体の挙動を解析するシミュレーション方法が知られている。このシミュレーション方法は、粒子法と呼ばれる。粒子法では、流体が複数の粒子で表され、複数の粒子のそれぞれについて、流体の支配方程式を数値的に解くことによって、流体の挙動を解析する。
【先行技術文献】
【非特許文献】
【0003】
【非特許文献1】原田隆宏、越塚誠一、「ＳＰＨにおける壁境界計算手法の改良」、情報処理学会論文誌、Ｖｏｌ．４８、Ｎｏ．４、１８３８～１８４６頁、２００７年
【発明の概要】
【発明が解決しようとする課題】
【0004】
流体機械の中には、空気、水、油等が混じった流れ、すなわち混相流を生じるものがある。混相流を数値的に解析するためには、各相の界面、例えば気液界面を適切に取り扱う必要がある。そのための手法は、従来から広く研究されており、レイノルズ数が比較的低い遅い流れにおいては、精度良く解析できることが分かっている。しかしながら、レイノルズ数が高い速い流れにおいては、界面が大きく変形したり飛沫を発生したりすることがあるため、与えられた空間解像度では界面を捉えることができなくなり、解析が不安定になることがある。
【0005】
図１Ａ及び図１Ｂを参照して、解析が不安定になる例について説明する。図１Ａ及び図１Ｂは、それぞれレイノルズ数が比較的低い気液混相流及びレイノルズ数が比較的高い気液混相流を表した複数の計算粒子１０の分布の一例を示す模式図である。液体を表す複数の計算粒子１０Ｌ及び気体を表す複数の計算粒子１０Ｇが混在している。図１Ａ及び図１Ｂにおいて、液体を表す計算粒子１０Ｌを灰色の円形で示し、気体を表す計算粒子１０Ｇを中空の円形で示している。
【0006】
レイノルズ数が比較的低い流れにおいては、図１Ａに示すように、気液界面の形状が単純である。これに対して、レイノルズ数が比較的高い流れにおいては、図１Ｂに示すように、気液界面の変形や飛沫の発生によって、界面の形状が複雑になる場合がある。界面の形状が、与えられた空間解像度で捉えきれなくなると、数値解析上では、局所的に物理量が不連続な領域が点在している状態になる。物理量が不連続な領域が存在すると、従来の方法では、支配方程式を解くために必要な場の物理量の導関数の値が極端に大きくなり、解析が数値的に不安定になってしまう。
【0007】
本発明の目的は、混相流体の挙動を安定に解析することが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。
【課題を解決するための手段】
【0008】
本発明の一観点によると、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーション方法であって、
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するシミュレーション方法が提供される。
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーション方法であって、
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するシミュレーション方法が提供される。
【0009】
本発明の他の観点によると、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーション方法であって、
前記運動量保存式は、解析における数値的な不安定性を低減させる人工粘性項を含んでおり、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるシミュレーション方法が提供される。
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーション方法であって、
前記運動量保存式は、解析における数値的な不安定性を低減させる人工粘性項を含んでおり、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるシミュレーション方法が提供される。
【0010】
本発明のさらに他の観点によると、
粒子法を用いて混相流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算する機能と
を有するシミュレーション装置が提供される。
粒子法を用いて混相流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算する機能と
を有するシミュレーション装置が提供される。
【0011】
本発明のさらに他の観点によると、
粒子法を用いて流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式、及び解析における数値的な不安定性を低減させる人工粘性項を含む運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす複数の計算粒子である参照粒子のそれぞれに付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える機能と
を有するシミュレーション装置が提供される。
粒子法を用いて流体の流れを解析するシミュレーション装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて演算を行う処理部と
を備え、
前記処理部は、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式、及び解析における数値的な不安定性を低減させる人工粘性項を含む運動量保存式を含む支配方程式を解いて時間発展させる機能と、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす複数の計算粒子である参照粒子のそれぞれに付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変える機能と
を有するシミュレーション装置が提供される。
【0012】
本発明のさらに他の観点によると、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーションの手順をコンピュータに実行させるプログラムであって、
前記複数の計算粒子のそれぞれについて質量保存式を解く手順において、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するプログラムが提供される。
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式及び運動量保存式を含む支配方程式を解いて時間発展させるシミュレーションの手順をコンピュータに実行させるプログラムであって、
前記複数の計算粒子のそれぞれについて質量保存式を解く手順において、計算対象の計算粒子である対象粒子に付与された密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記対象粒子に対する前記参照粒子の相対速度を小さくして計算するプログラムが提供される。
【0013】
本発明のさらに他の観点によると、
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式、及び解析における数値的な不安定性を低減させる人工粘性項を含む運動量保存式を含む支配方程式を解いて時間発展させるシミュレーションの手順をコンピュータに実行させるプログラムであって、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く手順において、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるプログラムが提供される。
混相流体を複数の計算粒子で表し、前記複数の計算粒子のそれぞれに物理量を付与し、解析空間に前記複数の計算粒子を配置し、前記複数の計算粒子のそれぞれに付与された物理量及び前記複数の計算粒子の位置を、質量保存式、及び解析における数値的な不安定性を低減させる人工粘性項を含む運動量保存式を含む支配方程式を解いて時間発展させるシミュレーションの手順をコンピュータに実行させるプログラムであって、
前記複数の計算粒子のそれぞれについて前記運動量保存式を解く手順において、計算対象の計算粒子である対象粒子の密度と、前記対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、前記人工粘性項による人工粘性の強さを変えるプログラムが提供される。
【発明の効果】
【0014】
複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、対象粒子に対する参照粒子の相対速度を小さくして計算することにより、密度が異なる計算粒子が混在する混相流体の挙動を安定に解析することが可能になる。
【0015】
また、複数の計算粒子のそれぞれについて運動量保存式を解く際に、計算対象の計算粒子である対象粒子の密度と、対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、運動量保存式の人工粘性項による人工粘性の強さを変えることにより、密度が異なる計算粒子が混在する混相流体の挙動を安定に解析することが可能になる。
【図面の簡単な説明】
【0016】
【発明を実施するための形態】
【0017】
本発明の一実施例について説明する前に、一般的なＳＰＨ（Smoothed Particle Hydrodynamics）法について簡単に説明する。ＳＰＨ法は、連続体を複数の計算粒子の集合を用いて離散化し、流れ場の時間発展を計算するラグランジュ型の数値計算手法である。ＳＰＨ法では、粒子１個の質量分布として、連続かつ微分可能なカーネル関数を用いるとともに、カーネル関数の微分係数の重ね合わせを空間微分とする離散化モデルを用いる。
【0018】
ＳＰＨ法では、計算粒子の質量ｍはカーネル関数Ｗによってその周囲に分配され、連続的に分布しており、密度場は、その重ね合わせとして与えられる。すなわち、ｉ番目の計算粒子の位置ｘｉにおける密度ρｉは、以下の式で表される。
本明細書において、下付きの添え字ｉ及びｊは、計算粒子のそれぞれに付与された通し番号を表し、下付きの添え字ｉ、ｊが付された物理量は、それぞれｉ番目及びｊ番目の計算粒子に付された物理量を意味する。
【数1】
【0019】
カーネル関数Ｗとして、例えば、以下の関数を用いることができる。以下の関数は、Ｃ２Ｗｅｎｄｌａｎｄカーネルと呼ばれる。
【数2】
【0020】
ｈは、平滑化長さであり、カーネル幅ともいわれる。カーネル幅ｈは、複数の計算粒子を同一の質量及び物性密度の球体に置き換えた等価球の直径に等しい。本明細書において、カーネル幅ｈを粒子の直径という場合がある。
【0021】
任意の位置ｘｉにおける任意のスカラー量Ａ及びベクトル量Ａの値は、カーネル関数の重ね合わせによって以下の式で表すことができる。
【数3】
【0022】
任意のスカラー量Ａの勾配は、カーネル関数の微分係数の重ね合わせによって以下の式で表される。
【数4】
【数5】
【0023】
ここで、∇Ｗ（ｘｉｊ，ｈ）はカーネル関数の導関数であり、カーネル関数としてＣ２Ｗｅｎｄｌａｎｄカーネルを用いる場合、以下の式で表される。
なお、式（２）の関数Ｗ（ｑ）は１変数の関数であるため、その勾配∇Ｗはスカラー量である。
【数6】
【0024】
任意のベクトル量Ａに対する発散は、式（４）と同様に、以下の式で表すことができる。
【数7】
【0025】
任意のスカラー量Ａに対する２次導関数は、例えば以下の公式を用いて算出することができる。
【数8】
【0026】
【0027】
図２は、本実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部３０、処理部３１、出力部３２、及び記憶部３３を含む。入力部３０から処理部３１にシミュレーション条件等が入力される。さらに、ユーザから入力部３０に各種指令（コマンド）等が入力される。入力部３０は、例えば通信装置、リムーバブルメディア読取装置、キーボード、ポインティングデバイス等で構成される。
【0028】
処理部３１は、入力されたシミュレーション条件及び指令に基づいて粒子法、例えばＳＰＨ法によるシミュレーションを実行する。さらに、シミュレーション結果を出力部３２に出力する。シミュレーション結果には、シミュレーション対象物である粒子系の粒子の状態、粒子系の物理量の時間的変化等を表す情報が含まれる。処理部３１は、例えばコンピュータの中央処理ユニット（ＣＰＵ）を含む。ＳＰＨ法によるシミュレーションをコンピュータに実行させるためのプログラムが、記憶部３３に記憶されている。出力部３２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。
【0029】
［流体層の支配方程式］
次に、本実施例によるシミュレーション方法で用いられる流体相の支配方程式について説明する。
次に、本実施例によるシミュレーション方法で用いられる流体相の支配方程式について説明する。
【0030】
流体相の支配方程式は、以下に示すナビエストークス（Navier-Stokes）方程式で記述される。
ここで、ｕは速度ベクトル、ｔは時間、ρは密度、ｐは圧力、及びμは粘性係数を表す。ｇは重力ベクトルである。
【数9】
【数10】
【0031】
支配方程式（９）、（１０）を、ＳＰＨ法を用いて離散化すると、例えば以下のように表すことができる。
ここで、τは粘性応力テンソルを表す。式（１４）は、流体が弱い圧縮性を持つという仮定により導出される状態方程式である。ｃは仮想的な音速、ρ０は基準密度を表す。音速ｃの値は、例えば流れ場の代表速度の１０倍程度に設定することができる。基準密度ρ０は、例えば作動流体の密度に設定することができる。式（１１）～式（１４）に示す支配方程式は、非圧縮性流れを、ＳＰＨ法を用いて解析する場合に広く用いられている。式（１２）は、質量保存式であり、式（１３）は、運動量保存式である。
【数11】
【数12】
【数13】
【数14】
【0032】
［混相流体の解析］
密度が大きく異なる２種類の流体が混在する流れ場を解析する際には、質量保存式（１２）を変形して、密度に関する式ではなく、体積に関する式にすると、解析上取り扱いが容易になる。式（１２）を、Ｖｉ＝ｍｉ／ρｉを用いて変形すると、以下の式が得られる。
密度が大きく異なる２種類の流体が混在する流れ場を解析する際には、質量保存式（１２）を変形して、密度に関する式ではなく、体積に関する式にすると、解析上取り扱いが容易になる。式（１２）を、Ｖｉ＝ｍｉ／ρｉを用いて変形すると、以下の式が得られる。
【数15】
【数16】
【0033】
【数17】
【数18】
【0034】
さらに、状態方程式（１４）は以下のように書き換えることができる。
ここで、Ｖ０は計算粒子の基準体積を表す。
【数19】
【0035】
［混相流体の解析における質量保存式の取り扱い］
次に、図３Ａ～図３Ｃを参照して、混相流体の解析における質量保存式の取り扱いについて説明する。
質量保存式（１８）の右辺の和を計算する際には、ｉ番目の計算粒子とｊ番目の計算粒子との密度差を考慮しなければ質量保存式を正しく満足することができない。以下、計算対象であるｉ番目の計算粒子を「対象粒子」といい、対象粒子に作用を及ぼす周囲のｊ番目の計算粒子を「参照粒子」ということとする。
次に、図３Ａ～図３Ｃを参照して、混相流体の解析における質量保存式の取り扱いについて説明する。
質量保存式（１８）の右辺の和を計算する際には、ｉ番目の計算粒子とｊ番目の計算粒子との密度差を考慮しなければ質量保存式を正しく満足することができない。以下、計算対象であるｉ番目の計算粒子を「対象粒子」といい、対象粒子に作用を及ぼす周囲のｊ番目の計算粒子を「参照粒子」ということとする。
【0036】
図３Ａは、解析空間に分布する気相（例えば空気）の物理量が付与された複数の計算粒子１０Ｇ（以下、単に「気相の計算粒子」という場合がある。）、及び液相（例えば水）の物理量が付与された複数の計算粒子１０Ｌ（以下、単に「液相の計算粒子」という場合がある。）の分布を示す模式図である。図３Ａにおいて、液相の計算粒子１０Ｌを灰色の円形で表し、気相の計算粒子１０Ｇを中空の円形で表す。図３Ｂ及び図３Ｃにおいても同様である。
【0037】
【0038】
図３Ｂは、対象粒子１０Ｇｉが気相の計算粒子である場合の一例を示す模式図である。図３Ｂにおいて、対象粒子１０Ｇｉを相対的に太い輪郭線で表している。対象粒子１０Ｇｉの周囲に分布する複数の気相の参照粒子１０Ｇｊ及び複数の液相の参照粒子１０Ｌｊについて、式（１８）の右辺の和を計算する。対象粒子１０Ｇｉには、速度ベクトルｕｉが付与されている。このとき、気相の参照粒子１０Ｇｊ及び液相の参照粒子１０Ｌｊのそれぞれに付与されている速度ベクトルｕｊをそのまま用いる。
【0039】
図３Ｃは、対象粒子１０Ｌｉが液相の計算粒子である場合の一例を示す模式図である。図３Ｃにおいて、対象粒子１０Ｌｉを相対的に太い輪郭線で表している。この場合、対象粒子１０Ｌｉへの気相の参照粒子１０Ｇｊの寄与を無視する。具体的には、気相の参照粒子１０Ｇｊの速度ベクトルｕｊの大きさ及び向きを、対象粒子１０Ｌｉの速度ベクトルｕｉの大きさ及び向きと等しくして、式（１８）の右辺の和を計算する。言い換えると、対象粒子１０Ｌｉに対する参照粒子１０Ｇｊの相対速度をゼロにする。
【0040】
【0041】
【数20】
【数21】
【数22】
【数23】
【0042】
例えば、気相の対象粒子１０Ｇｉと液相の参照粒子１０Ｌｊとについて計算するときは、ρｉ＜＜ρｊが成り立つため、式（２３）のｕｊ
＊－ｕｉは、ｕｊ－ｕｉとほぼ等しくなる。これは、参照粒子１０Ｇｊの速度ベクトルｕｊをそのまま用いて質量保存式を計算することを意味する。
【0043】
液相の対象粒子１０Ｌｉと気相の参照粒子１０Ｇｊとについて計算するときは、式（２３）のｕｊ
＊－ｕｉは、ｕｊ－ｕｉのρｊ／ρｉ倍になる。これは、対象粒子１０Ｌｉに対する参照粒子１０Ｇｊの相対速度ｕｊ－ｕｉを、対象粒子１０Ｌｉの密度ρｉに対する参照粒子１０Ｇｊの密度ρｊの比ρｊ／ρｉに応じて小さくすることを意味する。一般的に、ρｉ＞＞ρｊが成り立つため、相対速度ｕｊ－ｕｉはほぼゼロになる。
【0044】
ρｉ及びρｊの一方が他方より十分大きいという条件が満足されない場合は、式（２０）に示す密度の比に基づいて、対象粒子に対する参照粒子の相対速度ｕｊ－ｕｉを小さくするとよい。対象粒子の密度ρｉに対する参照粒子の密度ρｊの比ρｊ／ρｉが小さくなるほど、相対速度ｕｊ－ｕｉをより小さくするとよい。
【0045】
［混相流体の解析における運動量保存式の取り扱い］
次に、混相流体の解析における運動量保存式の取り扱いについて説明する。
密度差の大きな２つの流体の混合流れを解析する際には、２つの流体の界面において圧力の値が不連続になるため、解析が数値的に不安定になる場合がある。この不安定さを抑制するために、運動量保存式（１３）の右辺第１項（圧力勾配項）の和を計算する際に、圧力勾配が界面を隔てて不連続にならないようにする必要がある。そのための手法として、例えば圧力勾配項の和を計算する際に、着目する２個の計算粒子のそれぞれに付与された密度を用いて運動方程式（１３）の右辺第１項を変形する手法を用いることができる。
次に、混相流体の解析における運動量保存式の取り扱いについて説明する。
密度差の大きな２つの流体の混合流れを解析する際には、２つの流体の界面において圧力の値が不連続になるため、解析が数値的に不安定になる場合がある。この不安定さを抑制するために、運動量保存式（１３）の右辺第１項（圧力勾配項）の和を計算する際に、圧力勾配が界面を隔てて不連続にならないようにする必要がある。そのための手法として、例えば圧力勾配項の和を計算する際に、着目する２個の計算粒子のそれぞれに付与された密度を用いて運動方程式（１３）の右辺第１項を変形する手法を用いることができる。
【0046】
例えば、以下の式
を用いて、式（１３）の右辺第１項を以下のように変形する。
【数24】
【数25】
【0047】
運動方程式（１３）において、レイノルズ数の高い流れを解析する際には、右辺第１項の圧力に関する和の中に、解析における数値的な不安定性を低減させるために、人工粘性項を付加することができる。人工粘性項として、以下の式を用いることができる。
ここで、ζは人工粘性の強さを決めるパラメータであり、例えばζ＝０．１と設定することができる。κは計算粒子が重なっている場合においても分母の値がゼロにならないようにするための微小量であり、例えばκ＝０．０１と設定することができる。
【数26】
【0048】
レイノルズ数の高い混相流れにおいて、ζの値を大きくすることで、数値的に安定な解析が可能になる。ところが、ζを必要以上に大きくすると、流体が持つ物理的な粘性に対して人工粘性が卓越してしまい、流れ場が非物理的な挙動を示すことがある。対象とする流れにおいて数値的に不安定になりやすいのは２つの流体の界面近傍の計算粒子である。したがって、界面近傍で人工粘性が強く作用し、それ以外の領域では人工粘性が弱く作用するようにすれば、流れ場全体で物理的な粘性に対して人工粘性が卓越することがなくなる。
【0049】
界面近傍に存在する計算粒子とそうでない計算粒子に対して、人工粘性の強さを調節するために、人工粘性を以下の式で定義するとよい。
ここで、ａ１、ａ２は、それぞれ界面から離れた領域及び界面近傍の領域における人工粘性の強さを決めるパラメータであり、例えばａ１＝０．１、ａ２＝１．０と設定することができる。γ＊は、着目する計算粒子が界面近傍の領域に存在するか否かを判定するためのパラメータであり、例えばγ＊＝０．９５と設定することができる。
【数27】
【数28】
【0050】
ｂ１は、着目する計算粒子に付与されている密度に対して、それと異なる密度が付与された計算粒子が周囲にどの程度存在するかを評価するための指標であり、以下の式で定義することができる。
密度ρｉの対象粒子の周囲に存在するすべての参照粒子の密度ρｊが、対象粒子の密度ρｉと等しい場合、指標ｂ１＝１になる。対象粒子の周囲に、密度が異なる参照粒子の個数が増えるにしたがって、指標ｂ１は小さくなる。さらに、密度の比が１から離れるにしたがって、指標ｂ１は小さくなる。すなわち、密度の比が１から大きく離れている参照粒子が、周囲に多く存在するほど、指標ｂ１が小さくなる。密度の比が１から大きく離れている参照粒子が、周囲に多く存在するということは、対象粒子が、２つの流体の界面の近傍に位置することを意味する。
【数29】
【0051】
指標ｂ１がγ＊以下の場合に、対象粒子は２つの流体の界面の範囲内に位置すると判定される。このとき、式（２７）のγ＝ａ２＝１．０となり、人工粘性が強くなる。指標ｂ１がγ＊より大きい場合に、対象粒子は２つの流体の界面の範囲外に位置すると判定される。このとき、式（２７）のγ＝ａ１＝０．１となり、人工粘性が弱くなる。
【0052】
人工粘性に、式（２７）を用いることによって、密度の異なる計算粒子の分布に応じて、自動的に人工粘性の強さが調節される。
【0053】
式（２９）の計算結果は、対象粒子の密度ρｉと、複数の参照粒子のそれぞれの密度ρｊとの比に依存する。式（２９）の計算結果ｂ１によって、式（２８）で定義されるγの値が決まる。すなわち、式（２７）～式（２９）は、対象粒子に付与された密度と、複数の参照粒子のそれぞれに付与された密度との比に基づいて、人工粘性項による人工粘性の強さを変えることを意味する。
【0054】
［固相の支配方程式］
混相流体に粉粒体が混在する場合がある。次に、粉粒体の挙動を支配する固相の支配方程式について説明する。
混相流体に粉粒体が混在する場合がある。次に、粉粒体の挙動を支配する固相の支配方程式について説明する。
【0055】
固相の支配方程式は、以下のように記述することができる。なお、粉粒体を構成する複数の計算粒子に付される通し番号を下付き文字のａ、ｂで表し、流体を表す計算粒子に付与する番号ｉ、ｊと区別して表記する。
ここで、ｍ、Ｉ、ｖ、及びωは、それぞれ１個の計算粒子の質量、慣性モーメント、並進速度ベクトル、及び各速度ベクトルを表す。Ｆは力を表し、Ｔはトルクを表す。下付き文字ｆは、力またはトルクが流体から受けるものであることを意味し、下付き文字ｃは、力またはトルクが固相同士の相互作用によるものであることを意味する。
【数30】
【数31】
【数32】
【0056】
固相同士の相互作用による力Ｆｃ及びトルクＴｃの算出には、離散化要素法（ＤＥＭ法）で用いられる手法を適用することができる。また、粉粒体の計算粒子が流体から受けるトルクＴｆは、ここでは考慮しない。
【0057】
粉粒体の計算粒子が流体から受ける力は、圧力勾配力、粘性力、及び抗力の３つに分けることができ、以下の式で表すことができる。
ここで、τは粘性応力テンソル、Ｆｄは抗力を表す。圧力勾配力と粘性力とは、例えば以下の式を用いて算出することができる。
【数33】
【数34】
【0058】
抗力は、以下の式を用いて算出することができる。
ここで、εは、流体相の空隙率（１－固相の体積／流体相の体積）を表す。空隙率εの値は、流体相の計算粒子に付与し、次の式で算出する。
【数35】
【数36】
【0059】
粉粒体のａ番目の計算粒子の位置における空隙率は、以下の式で算出することができる。
【数37】
【0060】
式（３５）のｕａは、粉粒体のａ番目の粒子の位置における流体相の速度ベクトルを表し、例えば以下の式で算出することができる。
【数38】
【0061】
式（３５）のβａは以下の式で与えられる。
ここで、ρｆは、流体相の密度を表す。ＣＤは抗力係数であり、以下の式を用いて算出することができる。
ここで、Ｒｅａは、粒子直径を代表長さとするレイノルズ数であり、以下の式で表される。
【数39】
【数40】
【数41】
【0062】
流体相が固相から受ける力Ｓｉは、固相が流体相から受ける力（式（３３））の反作用として、以下の式を用いて算出することができる。
【数42】
【0063】
［計算粒子径未満の径を持つ泡のモデル］
水と油や、水と空気を攪拌した場合、無数の微小な泡が形成される。混相流体の流れを精度良く解析するためには、微細な泡の運動を再現することが好ましい。しかしながら、一般に泡の大きさは微小であるため、泡の形状を再現するためには極めて高い空間解像度が必要である。このため、微小な泡を含む大規模な流れを解析することは、計算資源の観点から困難である。限られた計算資源で泡の挙動を予測可能な解析手法が必要である。次に、図４を参照して、微小な泡の挙動を解析することが可能な泡モデルについて説明する。
水と油や、水と空気を攪拌した場合、無数の微小な泡が形成される。混相流体の流れを精度良く解析するためには、微細な泡の運動を再現することが好ましい。しかしながら、一般に泡の大きさは微小であるため、泡の形状を再現するためには極めて高い空間解像度が必要である。このため、微小な泡を含む大規模な流れを解析することは、計算資源の観点から困難である。限られた計算資源で泡の挙動を予測可能な解析手法が必要である。次に、図４を参照して、微小な泡の挙動を解析することが可能な泡モデルについて説明する。
【0064】
図４は、液相の物理量が付与された計算粒子１０Ｌの模式図である。液相の物理量が付与された計算粒子１０Ｌの中に、計算粒子１０Ｌの直径より小さな気相の複数の気泡１０ＬＢが存在していると考える。気泡１０ＬＢの直径をｄ、気相の体積分率をαと標記する。
【0065】
重力下においては、気相の粒子は浮力によって上昇する。このときの気相と液相との速度差をｕｒと標記すると、速度差ｕｒは以下の式で表現することができる。
ここで、νは液相の動粘性係数である。ＣＤは泡の抗力係数であり、以下のように書ける。
【数43】
【数44】
【0066】
Ｒｅは泡のレイノルズ数であり、以下の式で算出することができる。
ここで、ａＲ、ｂＲ、ｃＲは定数で、ａＲ＝０．６３０１２２６、ｂＲ＝２１２４５５２、ｃＲ＝２１０３７．８７である。また、Ｒρは、以下の式で計算できる。
ここで、ρｇａｓは気相の密度であり、ρｌｉｑは液相の密度である。
【数45】
【数46】
【数47】
【0067】
気相と液相との速度差ｕｒによって、液相に作用する力は以下のように表すことができる。
ここで、Ａは、１個の気泡１０ＬＢの投影面積であり、ｎは１個の計算粒子１０Ｌ内に存在する気泡の個数であり、Ｖｂは１個の気泡１０ＬＢの体積である。ｄは気泡１０ＬＢの直径であり、実験結果等に基づいて設定することができる。
【数48】
【数49】
【数50】
【数51】
【0068】
泡の移流によって、計算粒子１０Ｌのそれぞれの気相の体積分率が変化する。この現象を再現するためには、流体相の支配方程式に、以下の気相の体積分率に関する保存則の式を追加する必要がある。
ここで、Ｖαｉは、ｉ番目の計算粒子１０Ｌ内の気泡１０ＬＢの合計の体積であり、気相の体積分率αｉは、以下の式（５７）で表される。
【数52】
【0069】
［流体相の最終的な支配方程式］
上述の解析モデルをまとめると、流体相の最終的な支配方程式は以下のように表すことができる。
上述の解析モデルをまとめると、流体相の最終的な支配方程式は以下のように表すことができる。
【数53】
【数54】
【数55】
【数56】
【0070】
【数57】
【数58】
【数59】
【数60】
【0071】
上述の説明では、気相と液相との混合流体を考えたが、気相同士や液相同士の混合流体に対しても、同じ解析モデルを用いることができる。
【0072】
［シミュレーション方法のフローチャート］
次に、図５を参照して、本実施例によるシミュレーション方法について説明する。図５は、本実施例によるシミュレーション方法の手順を示すフローチャートである。各手順は、処理部３１（図２）によって実行される。
次に、図５を参照して、本実施例によるシミュレーション方法について説明する。図５は、本実施例によるシミュレーション方法の手順を示すフローチャートである。各手順は、処理部３１（図２）によって実行される。
【0073】
まず、処理部３１は、入力部３０に入力されたシミュレーション条件を取得する（ステップＳ１）。シミュレーション条件には、シミュレーション対象の流体を定義する情報、境界条件、及び初期条件等が含まれる。シミュレーション対象の流体を定義する情報には、混相流体を構成する２つの流体のそれぞれの密度、粘土等の物性値、及び粉粒体の物性値等が含まれる。境界条件には、解析空間の形状及び大きさを指定する情報、混相流体を表す計算粒子の境界近傍での取り扱い方法、及び粉粒体の境界近傍での取り扱い方法が含まれる。初期条件には、混相流体を表す計算粒子の初期分布及び初期物性値を指定する方法、粉粒体を表す計算粒子の初期分布及び初期物性値を指定する方法、及び複数の計算粒子を解析空間に配置するための情報が含まれる。
【0074】
処理部３１は、シミュレーション条件を取得すると、シミュレーション条件に基づいて、解析空間に混相流体及び粉粒体の複数の計算粒子を配置する（ステップＳ２）。
【0075】
次に、粉粒体を表す計算粒子の分布に基づいて、混相流体を表す計算粒子のそれぞれに対して、式（３６）を用いて空隙率を計算する（ステップＳ３）。次に、粉粒体を表す計算粒子のそれぞれについて、混相流体の計算粒子から受ける力を、式（３３）を用いて計算する（ステップＳ４）。次に、混相流体を表す計算粒子のそれぞれについて、粉粒体の計算粒子から受ける力を、式（４２）を用いて計算する（ステップＳ５）。
【0076】
混相流体の計算粒子のそれぞれについて、支配方程式（式（５３）～式（６０））を解くことによって、計算粒子のそれぞれに付与される物理量及び位置を更新する（ステップＳ６）。次に、粉粒体の計算粒子のそれぞれについて、支配方程式（式（３０）～式（３２））を解くことによって、計算粒子のそれぞれに付与される物理量及び位置を更新する（ステップＳ７）。
【0077】
ステップＳ３からステップＳ７までの手順を、計算終了条件が満たされるまで繰り返す（ステップＳ８）。計算終了条件が見たされると、処理部３１は出力部３２に計算結果を出力する（ステップＳ９）。
【0078】
次に、本実施例の優れた効果について説明する。
本実施例では、混相流体の複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、対象粒子に対する参照粒子の相対速度を小さくして計算する。これにより、解析において、２つの流体の密度差を考慮して質量保存則を満たすことができる。
本実施例では、混相流体の複数の計算粒子のそれぞれについて質量保存式を解く際に、計算対象の計算粒子である対象粒子に付与された密度と、対象粒子に作用を及ぼす他の複数の計算粒子である参照粒子に付与された密度との比に基づいて、対象粒子に対する参照粒子の相対速度を小さくして計算する。これにより、解析において、２つの流体の密度差を考慮して質量保存則を満たすことができる。
【0079】
また、本実施例では、複数の計算粒子のそれぞれについて、人工粘性項を含む運動量保存式を解く際に、対象粒子が、密度の異なる２つの流体の界面の範囲内に位置する場合、範囲外に位置する場合と比べて、人工粘性項による人工粘性の強さを大きくしている。これにより、界面近傍における数値解析の不安定性を抑制することができる。さらに、界面近傍以外の範囲では、物理的な粘性に対して人工粘性が卓越することによって流れ場が非物理的な挙動を示すことが抑制される。
【0080】
【0081】
図６は、解析モデルを示す模式図である。水と空気との混相流体、及び粉粒体の連成解析を行った。解析空間２０として、幅１００ｍｍ、高さ１００ｍｍ、奥行き２５ｍｍの直方体の空間を定義した。座標系として、幅方向をｘ軸、高さ方向をｙ軸、奥行き方向をｚ軸に設定した右手系のｘｙｚ直交座標系を用いる。初期状態において、解析空間２０の上半分の領域２０Ｕを、空気４０ｖｏｌ％、水６０ｖｏｌ％の混合流体で満たし、下半分の領域２０Ｌを、空気６０ｖｏｌ％、水４０ｖｏｌ％の混合流体で満たした。図６では、混相流体の計算粒子は示されていない。このとき、空気は、直径０．８４７ｍｍの泡として存在しているものとした。
【0082】
さらに、密度の異なる直径５ｍｍの粉粒体の計算粒子１０Ｓを、解析空間２０の奥行き方向の中心面に、高さ方向及び幅方向に１０ｍｍの間隔で等間隔に配置した。粉粒体の計算粒子１０Ｓの密度ρａを、初期状態における計算粒子１０Ｓの位置ベクトルのｘ方向成分ｘａを用いて、以下の式で定義した。
すなわち、計算粒子１０Ｓのそれぞれの密度は、初期位置においてｘ軸の正の向きに向かって徐々に高くなる。
【数61】
【0083】
重力加速度は９．８１ｍ／ｓ２とした。混合流体を表す計算粒子のカーネル幅は５ｍｍとした。時間の刻み幅を１．８７５×１０－５ｓとして、時刻０ｓから４ｓまで時間発展させた。
【0084】
図７ＡＡ～図８ＤＣは、解析によって得られた計算粒子の分布の時間変化を示す模式図である。図７ＡＡ～図７ＡＣは、初期状態（ｔ＝０）の分布を表し、図７ＢＡ～図７ＢＣは、０．５秒経過時点（ｔ＝０．５ｓ）の分布を表し、図８ＣＡ～図８ＣＣは、２秒経過時点（ｔ＝２ｓ）の分布を表し、図８ＤＡ～図８ＤＣは、４秒経過時点（ｔ＝４ｓ）の分布を表す。図７ＡＡ、図７ＢＡ、図８ＣＡ、及び図８ＤＡは、混相流体の計算粒子と粉粒体の計算粒子との分布を示している。図７ＡＢ、図７ＢＢ、図８ＣＢ、及び図８ＤＢは、混相流体の計算粒子のみの分布を示している。図７ＡＣ、図７ＢＣ、図８ＣＣ、及び図８ＤＣは、粉粒体の計算粒子のみの分布を示している。
【0085】
これらの模式図において、混相流体の計算粒子の濃度の違いは密度の違いを表しており、高密度の計算粒子を、低密度の計算粒子より相対的に濃く表している。図７ＡＢでは、上半分の計算粒子の密度が、下半分の計算粒子の密度より高いが、その差がわずかであるため、濃淡の差もわずかである。初期状態（図７ＡＢ）では、上半分に存在する流体の密度が下半分に存在する流体の密度より高いため、０．５秒経過した時点（図７ＢＢ）で、上半分に存在していた流体が密度差によって下方向に移動し始めていることがわかる。
【0086】
２秒経過した時点（図８ＣＢ）では、気泡として存在している空気が浮力によって上昇することにより、水と空気とが分離していく。解析空間の上部に存在する計算粒子の密度が空気の密度に近づき、下部に存在する計算粒子の密度が水の密度に近づく。４秒経過した時点（図８ＤＢ）では、水と空気とがほとんど完全に分離し、解析空間の上半分が空気、下半分が水で満たされていることがわかる。
【0087】
粉粒体の計算粒子に着目すると、０．５秒経過した時点（図７ＢＣ）では、水と空気とがほとんど分離していないため、水と空気との混合流体の密度より低い密度を持つ計算粒子が上昇し、それ以外の計算粒子が沈殿していることがわかる。２秒経過した時点（図８ＣＣ）では、水と空気とが分離するにつれて、０．５秒経過時点で沈殿していた計算粒子のうち、水と空気との混合流体の密度より高く、水の密度より低い密度を有する計算粒子が、水の表面の近傍に分布していることがわかる。また、０．５秒経過した時点で解析空間の上面の近傍に存在していた計算粒子が、空気と混合流体との界面の近傍まで下降していることがわかる。
【0088】
４秒経過した時点（図８ＤＣ）では、水の密度より高い密度の計算粒子が解析空間の下面に沈殿し、水の密度より低い密度の計算粒子が、水と空気との界面に分布していることがわかる。
【0089】
これらの解析結果から、本実施例によるシミュレーション方法を用いることにより、水と空気との混合流体と、粉粒体とが混在する流れ場の解析が可能であることが確認された。
【0090】
上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【0091】
１０ 計算粒子
１０Ｇ 気相の物理量が付与された計算粒子
１０Ｇｉ 気相の物理量が付与された対象粒子
１０Ｇｊ 気相の物理量が付与された参照粒子
１０Ｌｉ 液相の物理量が付与された対象粒子
１０Ｌｊ 液相の物理量が付与された参照粒子
１０Ｌ 液相の物理量が付与された計算粒子
１０ＬＢ 気泡
１０Ｓ 粉粒体の計算粒子
２０ 解析空間
２０Ｌ 解析空間の下半分の領域
２０Ｕ 解析空間の上半分の領域
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部
１０Ｇ 気相の物理量が付与された計算粒子
１０Ｇｉ 気相の物理量が付与された対象粒子
１０Ｇｊ 気相の物理量が付与された参照粒子
１０Ｌｉ 液相の物理量が付与された対象粒子
１０Ｌｊ 液相の物理量が付与された参照粒子
１０Ｌ 液相の物理量が付与された計算粒子
１０ＬＢ 気泡
１０Ｓ 粉粒体の計算粒子
２０ 解析空間
２０Ｌ 解析空間の下半分の領域
２０Ｕ 解析空間の上半分の領域
３０ 入力部
３１ 処理部
３２ 出力部
３３ 記憶部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
