(57)【要約】【課題】粒子の流動状態を解析する場合に、従来、個別要素法で再現していた現象を、連続体アプローチによっても再現しつつ、計算時間を短縮する。【解決手段】ＣＰＵ１１は、関連粘性係数として０でない定数を設定し、粒子径と、粒子間の摩擦特性を表す係数とを設定し、密度と、並行運動速度と、回転運動速度の初期値を設定して、それぞれＲＡＭ１３に記憶する。ＣＰＵ１１は、上記の係数と、密度と、並行運動速度と、回転運動速度と、並行運動速度および回転運動速度のそれぞれの空間微分値とを用いて、並行粘性係数を算出するとともに、密度と、この並行粘性係数を用いて、回転粘性係数を算出する。そして、ＣＰＵ１１は、設定された値や算出された値を用いて、密度の時間変化と、並行運動速度の時間変化と、回転運動速度の時間変化とをそれぞれ算出し、これら３つの状態量の所定時間後における値を算出する。【選択図】図１
2010003144.tif 000002
【特許請求の範囲】
【請求項１】
空間を複数に離散化した離散空間のそれぞれについて、当該離散空間に存在する粒子の質量の和を当該離散空間の体積で除算した値である密度と、当該離散空間に存在する各粒子が並行移動する速度の平均である並行運動速度と、当該離散空間に存在する各粒子の回転運動する速度の平均である回転運動速度と、前記粒子の径である粒子径と、前記粒子間の摩擦特性を表す係数と、前記並行運動速度と前記回転運動速度とが関連して生み出す応力に起因する係数である関連粘性係数として０でない定数とを記憶手段に記憶させる記憶制御手段と、
前記粒子間の摩擦特性を表す係数と、前記密度と、前記並行運動速度と、前記回転運動速度とを前記記憶手段から読み出し、読み出した内容と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値とを用いて、当該並行運動速度の空間微分値と前記離散空間にかかる応力とを結びつける係数である並行粘性係数を算出する第１算出手段と、
前記記憶手段から読み出した前記密度と、前記第１算出手段により算出された並行粘性係数を用いて、前記回転運動速度の空間微分値と前記離散空間にかかる偶応力とを結びつける係数である回転粘性係数を算出する第２算出手段と、
前記記憶手段から読み出した前記密度及び前記並行運動速度を用いて、当該密度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の密度を算出する第３算出手段と、
前記記憶手段から読み出した前記密度、前記並行運動速度及び前記回転運動速度と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値と、前記第１算出手段によって算出された前記並行粘性係数と、前記記憶手段から読み出した前記関連粘性係数とを用いて、前記並行運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の並行運動速度を算出する第４算出手段と、
前記記憶手段から読み出した前記密度、前記粒子径及び前記回転運動速度と、当該回転運動速度の空間微分値と、前記第２算出手段によって算出された前記回転粘性係数とを用いて、前記回転運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の回転運動速度を算出する第５算出手段と、
前記第３算出手段により算出された前記所定時間後の密度、前記第４算出手段により算出された前記所定時間後の並行運動速度、又は前記第５算出手段により算出された前記所定時間後の回転運動速度を、前記粒子の流動状態の解析結果として出力する出力手段と
を具備することを特徴とするデータ解析装置。
【請求項２】
空間を複数に離散化した離散空間のそれぞれについて、当該離散空間に存在する粒子の質量の和を当該離散空間の体積で除算した値である密度ρと、当該離散空間に存在する各粒子が並行移動する速度の平均である並行運動速度ｖiと、当該離散空間に存在する各粒子の回転運動する速度の平均である回転運動速度ωiと、前記粒子の径である粒子径ａと、前記粒子間の摩擦特性を表す係数φと、前記並行運動速度ｖiと前記回転運動速度ωiとが関連して生み出す応力σに起因する係数である関連粘性係数μ2として０でない定数とを記憶手段に記憶させる記憶制御手段と、
前記並行運動速度ｖiの空間微分値と前記離散空間にかかる応力σとを結びつける係数である並行粘性係数μ1を、数式（１）に基づいて算出する第１算出手段と、
【数１】
2010003144.tif 000013
前記回転運動速度ωiの空間微分値と前記離散空間にかかる偶応力γとを結びつける係数である回転粘性係数μ3を、数式（２）に基づいて算出する第２算出手段と
【数２】
2010003144.tif 000014
数式（３）に基づいて、前記密度ρの時間変化を算出し、算出した当該時間変化を用いて所定時間後の密度ρを算出する第３算出手段と、
【数３】
2010003144.tif 000015
数式（４）に基づいて、前記並行運動速度ｖiの時間変化を算出し、算出した当該時間変化を用いて所定時間後の並行運動速度ｖiを算出する第４算出手段と、
【数４】
2010003144.tif 000016
数式（５）に基づいて、前記回転運動速度ωiの時間変化を算出し、算出した当該時間変化を用いて所定時間後の回転運動速度ωiを算出する第５算出手段と、
【数５】
2010003144.tif 000017
前記第３算出手段により算出された前記所定時間後の密度ρ、前記第４算出手段により算出された前記所定時間後の並行運動速度ｖi、又は前記第５算出手段により算出された前記所定時間後の回転運動速度ωiを、前記粒子の流動状態の解析結果として出力する出力手段と
を具備することを特徴とするデータ解析装置。
【請求項３】
制御手段と記憶手段と出力手段とを備えたデータ解析装置に用いられるデータ解析方法であって、
前記制御手段が、空間を複数に離散化した離散空間のそれぞれについて、当該離散空間に存在する粒子の質量の和を当該離散空間の体積で除算した値である密度と、当該離散空間に存在する各粒子が並行移動する速度の平均である並行運動速度と、当該離散空間に存在する各粒子の回転運動する速度の平均である回転運動速度と、前記粒子の径である粒子径と、前記粒子間の摩擦特性を表す係数と、前記並行運動速度と前記回転運動速度とが関連して生み出す応力に起因する係数である関連粘性係数として０でない定数とを前記記憶手段に記憶させる記憶制御過程と、
前記制御手段が、前記粒子間の摩擦特性を表す係数と、前記密度と、前記並行運動速度と、前記回転運動速度とを前記記憶手段から読み出し、読み出した内容と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値とを用いて、当該並行運動速度の空間微分値と前記離散空間にかかる応力とを結びつける係数である並行粘性係数を算出する第１算出過程と、
前記制御手段が、前記記憶手段から読み出した前記密度と、前記第１算出過程により算出された並行粘性係数を用いて、前記回転運動速度の空間微分値と前記離散空間にかかる偶応力とを結びつける係数である回転粘性係数を算出する第２算出過程と、
前記制御手段が、前記記憶手段から読み出した前記密度及び前記並行運動速度を用いて、当該密度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の密度を算出する第３算出過程と、
前記制御手段が、前記記憶手段から読み出した前記密度、前記並行運動速度及び前記回転運動速度と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値と、前記第１算出手段によって算出された前記並行粘性係数と、前記記憶手段から読み出した前記関連粘性係数とを用いて、前記並行運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の並行運動速度を算出する第４算出過程と、
前記制御手段が、前記記憶手段から読み出した前記密度、前記粒子径及び前記回転運動速度と、当該回転運動速度の空間微分値と、前記第２算出手段によって算出された前記回転粘性係数とを用いて、前記回転運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の回転運動速度を算出する第５算出過程と、
前記制御手段が、前記第３算出過程により算出された前記所定時間後の密度、前記第４算出過程により算出された前記所定時間後の並行運動速度、又は前記第５算出過程により算出された前記所定時間後の回転運動速度を、前記粒子の流動状態の解析結果として、前記出力手段から出力する出力過程と
を具備することを特徴とするデータ解析方法。
【請求項４】
記憶手段と出力手段とを備えたコンピュータを、
空間を複数に離散化した離散空間のそれぞれについて、当該離散空間に存在する粒子の質量の和を当該離散空間の体積で除算した値である密度と、当該離散空間に存在する各粒子が並行移動する速度の平均である並行運動速度と、当該離散空間に存在する各粒子の回転運動する速度の平均である回転運動速度と、前記粒子の径である粒子径と、前記粒子間の摩擦特性を表す係数と、前記並行運動速度と前記回転運動速度とが関連して生み出す応力に起因する係数である関連粘性係数として０でない定数とを、前記記憶手段に記憶させる記憶制御手段と、
前記粒子間の摩擦特性を表す係数と、前記密度と、前記並行運動速度と、前記回転運動速度とを前記記憶手段から読み出し、読み出した内容と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値とを用いて、当該並行運動速度の空間微分値と前記離散空間にかかる応力とを結びつける係数である並行粘性係数を算出する第１算出手段と、
前記記憶手段から読み出した前記密度と、前記第１算出手段により算出された並行粘性係数を用いて、前記回転運動速度の空間微分値と前記離散空間にかかる偶応力とを結びつける係数である回転粘性係数を算出する第２算出手段と、
前記記憶手段から読み出した前記密度及び前記並行運動速度を用いて、当該密度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の密度を算出する第３算出手段と、
前記記憶手段から読み出した前記密度、前記並行運動速度及び前記回転運動速度と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値と、前記第１算出手段によって算出された前記並行粘性係数と、前記記憶手段から読み出した前記関連粘性係数とを用いて、前記並行運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の並行運動速度を算出する第４算出手段と、
前記記憶手段から読み出した前記密度、前記粒子径及び前記回転運動速度と、当該回転運動速度の空間微分値と、前記第２算出手段によって算出された前記回転粘性係数とを用いて、前記回転運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の回転運動速度を算出する第５算出手段と、
前記第３算出手段により算出された前記所定時間後の密度、前記第４算出手段により算出された前記所定時間後の並行運動速度、又は前記第５算出手段により算出された前記所定時間後の回転運動速度を、前記粒子の流動状態の解析結果として前記出力手段に出力させる出力制御手段と
して機能させるためのデータ解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子の流動状態を解析するデータ解析装置、データ解析方法及びデータ解析プログラムの技術に関する。
【背景技術】
【０００２】
計算機によって粒子の流動状態を解析する手法は、一般に、個別要素法（DEM: Discrete Element Method）によるものと連続体アプローチによるものの２つに大別される。非特許文献１には、個別要素法により粒子の流動状態を解析するモデルが開示されている。このような個別要素法を用いると、粒子個々の挙動を計算しなければならないため、一般に多大な時間を要する。特許文献１では、個別要素法により高速に大規模な計算を行うために個別要素法に特化した計算用ボードを用いる技術が開示されている。
【０００３】
一方、連続体アプローチによる解析では、粒子を流体のような連続体として扱う。連続体アプローチでは、連続体とみなした空間を離散化した各微小空間について、密度等の状態量を解析するので、粒子個々の挙動を計算する必要はなく、個別要素法に比べて高速に解析を行うことができる。非特許文献２には、粒子を極性流体モデルとして表現する手法が開示されている。また、非特許文献３には、粒子と同様、離散体である流氷の数値シミュレーション予測の手法が開示されている。
【非特許文献１】Cundall, P.A. & Strack, O.D.L.(1979). "A discrete numerical model for granular assemblies"  Geotechnique 29, No.1, 47−65
【特許文献１】特開平９−２５９１０５号公報
【非特許文献２】金谷健一（１９７９），「粉状体の流動の基礎理論（第１報，非圧縮性の流れ）」，日本機械学会論文集Ｂ編，４５巻３９２号，ｐ．５０７−５１４
【非特許文献３】Gutfraind, R. & Savage, S. B. (1997), "Smoothed Particle Hydrodynamics for the Simulation of Broken-Ice Fields: Mohr-Coulomb-Type Rheology and Frictional Boundary Conditions", Journal of Computational Physics 134, 203−215
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかしながら、非特許文献２の手法では、例えば個別要素法であれば再現できる非流動部の挙動（安息角の形成など）が再現できないという問題があった。また、非特許文献３の手法では、非流動部の挙動を再現できるものの、粒子の粒径や回転を表現できないという問題があった。
【０００５】
本発明の課題は、粒子の流動状態を解析する場合に、従来、個別要素法で再現していた現象を、連続体アプローチによっても再現しつつ、計算時間を短縮することである。
【課題を解決するための手段】
【０００６】
上述した課題を解決するため、本発明に係るデータ解析装置は、空間を複数に離散化した離散空間のそれぞれについて、当該離散空間に存在する粒子の質量の和を当該離散空間の体積で除算した値である密度と、当該離散空間に存在する各粒子が並行移動する速度の平均である並行運動速度と、当該離散空間に存在する各粒子の回転運動する速度の平均である回転運動速度と、前記粒子の径である粒子径と、前記粒子間の摩擦特性を表す係数と、前記並行運動速度と前記回転運動速度とが関連して生み出す応力に起因する係数である関連粘性係数として０でない定数とを記憶手段に記憶させる記憶制御手段と、前記粒子間の摩擦特性を表す係数と、前記密度と、前記並行運動速度と、前記回転運動速度とを前記記憶手段から読み出し、読み出した内容と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値とを用いて、当該並行運動速度の空間微分値と前記離散空間にかかる応力とを結びつける係数である並行粘性係数を算出する第１算出手段と、前記記憶手段から読み出した前記密度と、前記第１算出手段により算出された並行粘性係数を用いて、前記回転運動速度の空間微分値と前記離散空間にかかる偶応力とを結びつける係数である回転粘性係数を算出する第２算出手段と、前記記憶手段から読み出した前記密度及び前記並行運動速度を用いて、当該密度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の密度を算出する第３算出手段と、前記記憶手段から読み出した前記密度、前記並行運動速度及び前記回転運動速度と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値と、前記第１算出手段によって算出された前記並行粘性係数と、前記記憶手段から読み出した前記関連粘性係数とを用いて、前記並行運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の並行運動速度を算出する第４算出手段と、前記記憶手段から読み出した前記密度、前記粒子径及び前記回転運動速度と、当該回転運動速度の空間微分値と、前記第２算出手段によって算出された前記回転粘性係数とを用いて、前記回転運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の回転運動速度を算出する第５算出手段と、前記第３算出手段により算出された前記所定時間後の密度、前記第４算出手段により算出された前記所定時間後の並行運動速度、又は前記第５算出手段により算出された前記所定時間後の回転運動速度を、前記粒子の流動状態の解析結果として出力する出力手段とを具備することを特徴とする。
【０００７】
また、本発明に係るデータ解析方法は、制御手段と記憶手段と出力手段とを備えたデータ解析装置に用いられるデータ解析方法であって、前記制御手段が、空間を複数に離散化した離散空間のそれぞれについて、当該離散空間に存在する粒子の質量の和を当該離散空間の体積で除算した値である密度と、当該離散空間に存在する各粒子が並行移動する速度の平均である並行運動速度と、当該離散空間に存在する各粒子の回転運動する速度の平均である回転運動速度と、前記粒子の径である粒子径と、前記粒子間の摩擦特性を表す係数と、前記並行運動速度と前記回転運動速度とが関連して生み出す応力に起因する係数である関連粘性係数として０でない定数とを前記記憶手段に記憶させる記憶制御過程と、前記制御手段が、前記粒子間の摩擦特性を表す係数と、前記密度と、前記並行運動速度と、前記回転運動速度とを前記記憶手段から読み出し、読み出した内容と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値とを用いて、当該並行運動速度の空間微分値と前記離散空間にかかる応力とを結びつける係数である並行粘性係数を算出する第１算出過程と、前記制御手段が、前記記憶手段から読み出した前記密度と、前記第１算出過程により算出された並行粘性係数を用いて、前記回転運動速度の空間微分値と前記離散空間にかかる偶応力とを結びつける係数である回転粘性係数を算出する第２算出過程と、前記制御手段が、前記記憶手段から読み出した前記密度及び前記並行運動速度を用いて、当該密度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の密度を算出する第３算出過程と、前記制御手段が、前記記憶手段から読み出した前記密度、前記並行運動速度及び前記回転運動速度と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値と、前記第１算出手段によって算出された前記並行粘性係数と、前記記憶手段から読み出した前記関連粘性係数とを用いて、前記並行運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の並行運動速度を算出する第４算出過程と、前記制御手段が、前記記憶手段から読み出した前記密度、前記粒子径及び前記回転運動速度と、当該回転運動速度の空間微分値と、前記第２算出手段によって算出された前記回転粘性係数とを用いて、前記回転運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の回転運動速度を算出する第５算出過程と、前記制御手段が、前記第３算出過程により算出された前記所定時間後の密度、前記第４算出過程により算出された前記所定時間後の並行運動速度、又は前記第５算出過程により算出された前記所定時間後の回転運動速度を、前記粒子の流動状態の解析結果として、前記出力手段から出力する出力過程とを具備することを特徴とする。
【０００８】
また、本発明に係るプログラムは、記憶手段と出力手段とを備えたコンピュータを、空間を複数に離散化した離散空間のそれぞれについて、当該離散空間に存在する粒子の質量の和を当該離散空間の体積で除算した値である密度と、当該離散空間に存在する各粒子が並行移動する速度の平均である並行運動速度と、当該離散空間に存在する各粒子の回転運動する速度の平均である回転運動速度と、前記粒子の径である粒子径と、前記粒子間の摩擦特性を表す係数と、前記並行運動速度と前記回転運動速度とが関連して生み出す応力に起因する係数である関連粘性係数として０でない定数とを、前記記憶手段に記憶させる記憶制御手段と、前記粒子間の摩擦特性を表す係数と、前記密度と、前記並行運動速度と、前記回転運動速度とを前記記憶手段から読み出し、読み出した内容と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値とを用いて、当該並行運動速度の空間微分値と前記離散空間にかかる応力とを結びつける係数である並行粘性係数を算出する第１算出手段と、前記記憶手段から読み出した前記密度と、前記第１算出手段により算出された並行粘性係数を用いて、前記回転運動速度の空間微分値と前記離散空間にかかる偶応力とを結びつける係数である回転粘性係数を算出する第２算出手段と、前記記憶手段から読み出した前記密度及び前記並行運動速度を用いて、当該密度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の密度を算出する第３算出手段と、前記記憶手段から読み出した前記密度、前記並行運動速度及び前記回転運動速度と、当該並行運動速度及び当該回転運動速度のそれぞれの空間微分値と、前記第１算出手段によって算出された前記並行粘性係数と、前記記憶手段から読み出した前記関連粘性係数とを用いて、前記並行運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の並行運動速度を算出する第４算出手段と、前記記憶手段から読み出した前記密度、前記粒子径及び前記回転運動速度と、当該回転運動速度の空間微分値と、前記第２算出手段によって算出された前記回転粘性係数とを用いて、前記回転運動速度の時間変化を算出し、算出した当該時間変化を用いて所定時間後の回転運動速度を算出する第５算出手段と、前記第３算出手段により算出された前記所定時間後の密度、前記第４算出手段により算出された前記所定時間後の並行運動速度、又は前記第５算出手段により算出された前記所定時間後の回転運動速度を、前記粒子の流動状態の解析結果として前記出力手段に出力させる出力制御手段として機能させるためのデータ解析プログラムである。
【発明の効果】
【０００９】
本発明によれば、粒子の流動状態を解析するデータ解析装置において、従来、個別要素法で再現していた現象を、連続体アプローチによっても再現しつつ、計算時間を短縮することができる。
【発明を実施するための最良の形態】
【００１０】
以下、図面を参照して本発明の実施の形態を説明する。
（Ａ：構成）
（Ａ−１：装置構成）
図１は、実施形態に係るデータ解析装置１の全体構成を示すブロック図である。図１に示すように、データ解析装置１は、ＣＰＵ（Central Processing Unit）１１、ＲＯＭ（Read Only Memory）１２、ＲＡＭ（Random Access Memory）１３、記憶部１４、操作部１５及び表示部１６を備えている。これらは、バス１９に接続されている。
【００１１】
ＣＰＵ１１は、ＲＯＭ１２に記憶されているブートローダや記憶部１４に記憶されているコンピュータプログラムを読み出して実行することによりデータ解析装置１の各部を制御する制御手段である。ＲＯＭ１２は半導体素子等で構成された読み出し専用の不揮発性記憶装置である。このＲＯＭ１２には、記憶部１４に記憶されている基本ＯＳ（Operating System）をＲＡＭ１３にロードするためのブートローダやＢＩＯＳ（Basic Input/Output System）などが記憶されている。ＲＡＭ１３は半導体素子等で構成された揮発性の記憶手段であり、ＣＰＵ１１がプログラムを実行する際のワークエリアとして利用される。記憶部１４はハードディスクドライブなどの大容量の記憶手段であり、ＣＰＵ１１に読み込まれるプログラムを記憶する。
【００１２】
操作部１５は、キーボードやマウスなどを備えており、ユーザによる操作を受け付けてその操作内容に応じた信号（以下、操作信号という）をＣＰＵ１１に供給する。表示部１６は、ＣＰＵ１１の制御の下、種々の情報（文字、画像等）を表示する出力手段であり、例えば液晶ディスプレイである。
【００１３】
（Ａ−２：極性流体モデル）
次に、本実施形態においてデータ解析装置１が行うデータ解析の極性流体モデルを説明する。
【００１４】
粒子の挙動を解析する極性流体モデルでは、同一の剛体球である粒子が複数存在する空間を考える。この全ての粒子について、直径（以下、粒子径）はａであり、質量はｍである。この空間には重力が作用するが、粒子以外に質量を有するものはない。粒子同士は互いに接触摩擦による作用を及ぼしあうが、吸着などその他の作用を生じない。
【００１５】
この空間におけるある瞬間の状態量として、密度ρと並行運動速度ｖiと回転運動速度ωiとがある。密度ρとは、空間に存在する粒子の質量の和をその空間の体積で除算したものである。例えば、体積がＶの空間にｎ個の粒子が存在しているのであれば、その空間の密度ρは、ｎ×ｍ／Ｖで現される。並行運動速度ｖiとは、空間に存在する各粒子が並行移動する速度の平均（成分ごとの相加平均、以下同じ）であり、ベクトルで現される。回転運動速度ωiとは、空間に存在する各粒子が回転運動する速度の平均であり、ベクトルで現される。なお、以下、数式中においては、必要に応じてベクトル及びテンソルを添え字による表記と併用してボールド字体で表す。
【００１６】
データ解析装置１は、並行粘性係数μ1と、関連粘性係数μ2及び回転粘性係数μ3とを用いる。並行粘性係数μ1とは、ある空間における並行運動速度ｖiの空間微分値（勾配）とその空間にかかる応力とを結びつける係数である。並行粘性係数μ1は、次の数式（１）で現される。
【００１７】
【数１】
2010003144.tif 000003
【００１８】
ここで、Ｐは巨視的な圧力を示している。また、φは粒子間の摩擦特性を表す係数であり、ここでは安息角を示している。安息角とは、積み上げた粉体が自発的に崩れることなく安定を保つときの斜面の角度をいう。
また、ｆは極性流体の相当ひずみ速度を表している。このｆは並行運動速度ｖi、回転運動速度ωi及び粒子径ａを変数として定まる関数として現される。
【００１９】
関連粘性係数μ2とは、並行運動速度ｖiと回転運動速度ωiとが関連して生み出す応力に起因する係数である。
上述の非特許文献２では、関連粘性係数μ2が並行粘性係数μ1に比例すると定義されている。そのため、粒子の回転によるエネルギの散逸を計算で表すことができるが、非流動部の挙動（安息角の形成など）は再現できない。また、上述の非特許文献３では、関連粘性係数μ2を０としていたため、非流動部の挙動（安息角の形成など）は再現できるものの、粒子の回転によるエネルギの散逸を計算で表すことができない。
本実施形態に係るデータ解析装置１が用いる極性流体モデルでは、関連粘性係数μ2を並行粘性係数μ1とは独立の０でない定数としている。
【００２０】
回転粘性係数μ3とは、ある空間における回転運動速度ωiの勾配とその空間にかかる偶応力とを結びつける係数であり、Ｓｈｌｉｏｍｉｓの理論に従って、粒子の平均粒径ａと並行粘性係数μ1とから、次の数式（２）で現される。
【００２１】
【数２】
2010003144.tif 000004
【００２２】
密度ρ、並行運動速度ｖi、回転運動速度ωiのそれぞれの時間微分には、以下の数式（３）〜（５）が成り立つ。
【００２３】
【数３】
2010003144.tif 000005
【００２４】
【数４】
2010003144.tif 000006
【００２５】
【数５】
2010003144.tif 000007
【００２６】
数式（３）は、質量保存則に基づく質量保存方程式であり、ある空間の並行運動速度ｖiの発散に密度ρを乗算した値と、その空間の密度ρの時間微分との和が常に０であることを現している。
【００２７】
数式（４）は、運動量保存則に基づく運動方程式であり、空間にかかる単位体積当りの応力σの勾配が、加速度のベクトル（すなわち、並行運動速度ｖiの時間微分）と密度ρ（すなわち、単位体積当りの質量）との積に等しいことを示している。
この応力σは、密度ρ、並行運動速度ｖi、回転運動速度ωi、並行粘性係数μ1及び関連粘性係数μ2を変数として定まる関数σ（ρ，ｖi，ωi，μ1，μ2）として定義される。
【００２８】
数式（５）は、角運動量保存則に基づく角運動量方程式である。ここで、εijkは交代テンソルであり、γは偶応力を表している。偶応力γとは、単位面積あたりに作用する物体を回転させようとする力（トルク）を表す値であり、回転粘性係数μ3と回転運動速度ωiとを変数として定まる関数として現される。具体的には以下の数式（６）により現される。すなわち、偶応力γは回転粘性係数μ3と回転運動速度ωiの空間微分との積で現される。
【００２９】
【数６】
2010003144.tif 000008
【００３０】
以上により、本実施形態で用いる極性流体モデルでは、空間における密度ρ、並行運動速度ｖi及び回転運動速度ωiの時間微分を表す常微分方程式を、上述した数式（３）〜（５）として連立する。したがって適宜、時間と空間を離散化して、この連立常微分方程式に、密度ρ、並行運動速度ｖi及び回転運動速度ωiの初期値を設定し、ルンゲ・クッタ法や修正オイラー法などの数値積分法で積分することにより、これらの状態量の挙動・経時変化を解析することができる。
【００３１】
2010003144.tif 000009
【００３２】
【数７】
2010003144.tif 000010
2010003144.tif 000011
【００３３】
【数８】
2010003144.tif 000012
【００３４】
（Ｂ：動作）
次に、データ解析装置１の動作を説明する。
図２は、データ解析装置１の処理の流れを示すフロー図である。
データ解析装置１のＣＰＵ１１は、まず、後述する数値積分で用いる微小時間Δｔを設定するとともに（ステップＳＡ１０１）、関連粘性係数μ2を設定する（ステップＳＡ１０２）。そして、ＣＰＵ１１は、圧力Ｐ、安息角φ、平均粒子径ａを設定し（ステップＳＡ１０３）、計算対象である離散化された空間のそれぞれに対し、密度ρ、並行運動速度ｖi及び回転運動速度ωiの初期値を設定する（ステップＳＡ１０４）。これらの設定はＣＰＵ１１が、ＲＡＭ１３に記憶することにより行う。なお、ここで、ステップＳＡ１０１〜ＳＡ１０４は、この順番で実行しなくてもよい。また、上述した各設定値は、記憶部１４から読み出されてもよいし、操作部１５を介してユーザが指示した内容により決定されてもよい。
【００３５】
ＣＰＵ１１は、時刻ｔ（単位：秒）に初期値「０」を設定すると（ステップＳＡ１０５）、圧力Ｐ、安息角φ、並行運動速度ｖi、回転運動速度ωi及び平均粒子径ａをＲＡＭ１３から読み出し、読み出した内容と、並行運動速度ｖi及び回転運動速度ωiのそれぞれの空間微分値とを用いて、上述した数式（１）にしたがって並行粘性係数μ1を算出し、ＲＡＭ１３に記憶させる（ステップＳＡ１０６）。続いて、ＣＰＵ１１は、平均粒子径ａと、ステップＳＡ１０６で算出した並行粘性係数μ1を、これらを用いて、上述した数式（２）に従って、回転粘性係数μ3を算出し、ＲＡＭ１３に記憶させる（ステップＳＡ１０７）。次に、ＣＰＵ１１は、ＲＡＭ１３から読み出した密度ρ、並行運動速度ｖi、回転運動速度ωi、並行粘性係数μ1及び関連粘性係数μ2を用いて、応力σを算出するとともに（ステップＳＡ１０８）、回転粘性係数μ3及び回転運動速度ωiを用いて偶応力γを算出して（ステップＳＡ１０９）、これらをＲＡＭ１３に記憶させる。なお、ステップＳＡ１０７とＳＡ１０８は逆順に行ってもよく、ステップＳＡ１０８とＳＡ１０９は逆順に行ってもよい。要するに、ステップＳＡ１０７とＳＡ１０８に先行してステップＳＡ１０６が行われていればよく、ステップＳＡ１０９に先行してステップＳＡ１０７が行われていればよい。
【００３６】
並行粘性係数μ1、回転粘性係数μ3、応力σ及び偶応力γが算出され、ＲＡＭ１３に記憶された後、ＣＰＵ１１は、上述した数式（３）〜（５）を連立常微分方程式として、時刻ｔが微小時間Δｔだけ未来に進んだときの密度ρ、並行運動速度ｖi及び回転運動速度ωiを推算してこれらを更新する（ステップＳＡ１１０）。具体的にはＣＰＵ１１は、時刻ｔが微小時間Δｔだけ未来に進んだときの密度ρ、並行運動速度ｖi及び回転運動速度ωiをＲＡＭ１３に記憶するとともに、これらの全部又は一部を粒子の流動状態の解析結果として、表示部１６によって表示（出力）する。ＣＰＵ１１は、この推算に、例えばルンゲ・クッタ法や修正オイラー法などの数値積分法を用いる。そして、ＣＰＵ１１は、所定の終了条件が成立しているか否かを判断し（ステップＳＡ１１１）、終了条件が成立している場合には（ステップＳＡ１１１；ＹＥＳ）、データ解析を終了する。一方、終了条件が成立していない場合には（ステップＳＡ１１１；ＮＯ）、時刻ｔに微小時間Δｔを加算した値を、新たな時刻ｔとしてＲＡＭ１３に記憶し、処理をステップＳＡ１０６に戻す。このようにして、データ解析装置１は、所定の終了条件が満たされるまで、数値積分を行い、粒子の挙動・経時変化が計算される。
【００３７】
以上の通り、ＣＰＵ１１は、状態量の推算に際して、時刻ｔに対応する並行粘性係数μ1を逐次、算出するが、関連粘性係数μ2は一貫してステップＳＡ１０２で設定された定数を用いる。このため、シミュレーション全体を通じて、関連粘性係数μ2は並行粘性係数μ1に影響されることが無い。
【００３８】
そして、このデータ解析装置１は、例えば粉粒体流の一種である土石流や雪崩の流動状態について安息角を含めて計算できるので、これらを防ぐための防護柵を設計する際に有用である。また、複写機内のトナーの挙動解析を行うことでトナーが堆積する可能性が高い箇所を推定できるので、複写機の設計に役立てることができる。
【００３９】
（Ｃ：解析結果の例）
次に、本実施形態のデータ解析装置１を用いて、粒子の挙動をデータ解析した結果を、図３を用いて説明する。図３（ａ）、（ｂ）、（ｃ）、（ｄ）、（ｅ）及び（ｆ）は、シミュレーションにおいて、それぞれ時刻ｔ＝０，０．１３４６，０．２７２７，０．４０９１，０．５４５４，０．６８１８における粒子の挙動を示す図である。このシミュレーションにおいては、重力の方向に垂直な平面と、この平面から重力の反対方向に向かって延びる２枚の平板とによって囲われる空間を、２次元でモデル化している。そして、上述した平板の間に並行な壁を設けて、いずれか一方の平板とこの壁と平面で囲われる空間に粒子を充填させた場合を想定し、時刻ｔ＝０にこの壁を取り除くことによって起きる粒子の挙動を推算する。
【００４０】
図３（ａ）〜（ｆ）には、それぞれ上下２段の図を表している。上段の図は、各粒子の並行移動をベクトルとして表すとともに、その速度（０〜０．９９ｍ／ｓ）を青から赤まで段階的に区切られた色彩で表現している。下段の図は、各粒子の位置を点描によって表すとともに、その回転速度（−３．１４〜３．１４ｓ-1）を青から赤まで段階的に区切られた色彩で表現している。なお、ここでは図における反時計回りを回転速度の正方向として定義する。
【００４１】
図３（ａ）に示すように、初期状態である時刻ｔ＝０においては、粒子は壁と平板とに囲まれて静止している。壁が取り除かれて、時刻ｔ＝０．１３６４に至ると、図３（ｂ）に示すように、空間の左下隅に三角形の粒子群の層が現れる。この層は、図３（ｂ）の上段に示すように、並行移動をしておらず、下段に示すように、回転もしていない。ここでは、この層を不動層という。この不動層の上には、不動層の斜面に沿って並行移動する粒子群の層が形成されている。以下、この層を流動層という。流動層においては、図３（ｂ）〜（ｄ）の上段に示すように、不動層から離れるに従って、並行移動速度が上昇している。一方、回転速度については、図３（ｂ）〜（ｄ）の下段に示すように、流動層の内部の粒子群は時計回りにほぼ同じ速度で回転しつつ、粒子のない上側の空間との境界や、不動層との境界においては粒子群が、上述した内部よりも遅い速度で時計回りに回転していることが分かる。そして、図３（ｅ）（ｆ）に示すように、時刻ｔ＝０．５４５４以降は、粒子群は三角形の形状になり、略全体が静止状態に至ることが分かる。
なお、この図３のカラー画像を別途、手続補足書により提出する。
【００４２】
以上のとおり、本実施形態のデータ解析装置１では、関連粘性係数μ2に０でない定数を用いるため、非流動部の挙動を再現できるとともに、粒子の粒径や回転を表現することができる。また、本実施形態のデータ解析装置１は、連続体アプローチによってデータを解析しており、個別要素法のように個々の粒子の挙動を解析する必要が無いため、個別要素法を用いたデータ解析に比べて計算時間を短縮することができる。また、個別要素法を用いたデータ解析に比べると、処理負荷を抑えることができるので、専用の計算機を用いる必要はなく、汎用の計算機によってデータ解析を行うことができる。
【００４３】
（Ｄ：変形例）
上述の実施形態を以下のように変形してもよい。
（１）上述の実施形態においては、データ解析に、予め指定した安息角φを用いていたが、粒子間の摩擦特性を表す係数であれば、他の係数を用いてもよい。例えば、粉子の圧縮降伏応力や摩擦係数などを測定した値を用いてもよい。粉子体の圧縮降伏応力とは具体的には「粉体層に加重を加えた際に，固体状に静止している粉体層が流体状に流動を開始するときの粉体の内部応力値」である。また、粉子の摩擦係数とは具体的には「計測対象の粉体層に垂直方向力を与えてその粉体層を水平方向に移動させるのに必要な水平方向力とさきに加えた垂直力の比（すなわち，水平方向力／鉛直方向力）」である。なお、これら「粒子間の摩擦特性を表す係数」の具体的な計測方法としては，特開２００５−２４１４２５号公報や特開２００６−２６６７７０号公報に記載された計測方法がある。
【００４４】
（２）上述の実施形態においては、この全ての粒子について、粒子径と質量を一定値としたが、これらの値についてばらつきを考慮してもよい。この場合、粒子径や質量のばらつきの取り扱い方法として、平均値・中央値・最大値などをその粉粒体の代用特性として用いてもよいし、粒子径分布や粒子質量の分布に合わせて，個々の離散空間の粒子径又は質量を個々に定めてもよい。
【図面の簡単な説明】
【００４５】
【図１】実施形態に係るデータ解析装置１の全体構成を示すブロック図である。
【図２】データ解析装置１の処理の流れを示すフロー図である。
【図３（ａ）】時刻ｔ＝０における粒子の挙動の推算結果を示す図である。
【図３（ｂ）】時刻ｔ＝０．１３４６における粒子の挙動の推算結果を示す図である。
【図３（ｃ）】時刻ｔ＝０．２７２７における粒子の挙動の推算結果を示す図である。
【図３（ｄ）】時刻ｔ＝０．４０９１における粒子の挙動の推算結果を示す図である。
【図３（ｅ）】時刻ｔ＝０．５４５４における粒子の挙動の推算結果を示す図である。
【図３（ｆ）】時刻ｔ＝０．６８１８における粒子の挙動の推算結果を示す図である。
【符号の説明】
【００４６】
１…データ解析装置、１１…ＣＰＵ、１２…ＲＯＭ、１３…ＲＡＭ、１４…記憶部、１５…操作部、１６…表示部、１９…バス。
【図１】
2010003144.tif 000018
【図２】
2010003144.tif 000019
【図３（ａ）】
2010003144.tif 000020
【図３（ｂ）】
2010003144.tif 000021
【図３（ｃ）】
2010003144.tif 000022
【図３（ｄ）】
2010003144.tif 000023
【図３（ｅ）】
2010003144.tif 000024
【図３（ｆ）】
2010003144.tif 000025