シミュレーション方法および解析装置
【課題】流れのシミュレーションにおいて、より容易に粒子の動きを把握する。
【解決手段】解析装置１００は、複数の粒子を含む系を仮想空間内に定義する粒子定義部１０８と、粒子定義部１０８によって定義された系に、流れを表すよう制約を課す制約付与部１１０と、制約付与部１１０によって制約が課された系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部１２０と、を備える。数値演算部１２０は流れの擾乱源と関連付けられた粒子を特定する。本解析装置１００はさらに、特定された粒子を、そうでない粒子とは異なる態様でディスプレイ１０４に表示させる表示制御部１１８を備える。
【解決手段】解析装置１００は、複数の粒子を含む系を仮想空間内に定義する粒子定義部１０８と、粒子定義部１０８によって定義された系に、流れを表すよう制約を課す制約付与部１１０と、制約付与部１１０によって制約が課された系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部１２０と、を備える。数値演算部１２０は流れの擾乱源と関連付けられた粒子を特定する。本解析装置１００はさらに、特定された粒子を、そうでない粒子とは異なる態様でディスプレイ１０４に表示させる表示制御部１１８を備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子系を使用したシミュレーション方法および粒子系を解析する解析装置に関する。
【背景技術】
【０００２】
従来、古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）に基づくシミュレーションが知られている。ＭＤ法は物性をポテンシャルエネルギ関数で与えるのでモデル化の少ない手法である。
【０００３】
非特許文献１、２は、ＭＤ法による流体解析について報告している。非特許文献１、２では、ＭＤ法を用いてシリンダ周りの流れの計算が行われ、速度ベクトルの分布によってカルマン渦の可視化が試みられている。より具体的には、領域が分割され、分割された領域ごとの平均速度の向きが矢印で示されている。
【０００４】
また、ＭＤ法以外の粒子法（例えば、特許文献１、非特許文献３参照）では、速度分布以外にも、圧力分布や密度分布などによって渦を可視化することが試みられている。
有限要素法でも、速度分布や圧力分布以外にも渦度の分布から渦の確認が試みられている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００６−２８５８６６号公報
【非特許文献】
【０００６】
【非特許文献１】Ｄ．Ｃ．Ｒａｐａｐｏｒｔ，Ｅ．Ｃｌｅｍｅｎｔｉ、Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．、１９８６、ｖｏｌ５７、ｎｏ６、ｐ．６９５−６９８
【非特許文献２】Ｄ．Ｃ．Ｒａｐａｐｏｒｔ、Ｐｈｙｓ．Ｒｅｖ．Ａ、１９８７、ｖｏｌ３６、ｎｏ７、ｐ．３２８８−３２９９
【非特許文献３】越塚誠一、粒子法（計算力学レクチャーシリーズ５）、日本計算工学会編、丸善、２００５年
【非特許文献４】池畑光尚、船舶海洋工学のための流体力学入門、船舶技術協会、１９９３年、ｐ．１３０
【非特許文献５】白倉昌明、大橋秀雄、標準機械工学講座２３ 流体力学（２）、コロナ社、１９６９年、第２章
【非特許文献６】日本機械学会編、流れ−写真集、丸善、１９８４年、ｐ．４
【非特許文献７】丸山茂夫、井上和洋、「第１３回計算力学講演会講演論文集」、２０００年、ｐ．３８９−３９０
【非特許文献８】Ｄ．Ｌｅｖｅｓｑｕｅ，Ｌ．Ｖｅｒｌｅｔ，Ｊ．Ｋｕｒｋｉｊａｒｖｉ、Ｐｈｙｓ．Ｒｅｖ．Ａ、１９７３、ｖｏｌ７、ｎｏ５、ｐ．１６９０−１７００
【発明の概要】
【発明が解決しようとする課題】
【０００７】
上記の従来手法はいずれも、仮想空間内の単位領域ごとに速度や圧力などの物理量を計算し、計算された物理量をグラデーション等により見やすくしてユーザに提供することにより渦の可視化を目指している。このような手法は定常状態の解析には向いているかもしれないが、渦のように生来的に動的な現象をシミュレートする場合、流体粒子の動きを把握することが困難であり不向きである。
【０００８】
このような課題は渦のシミュレーションに限らず一般的な流れの擾乱のシミュレーションでも起こりうる。
【０００９】
本発明はこうした課題に鑑みてなされたものであり、その目的は、流れのシミュレーションにおいて、より容易に粒子の動きを把握することができる解析技術の提供にある。
【課題を解決するための手段】
【００１０】
本発明のある態様は、シミュレーション方法に関する。このシミュレーション方法は、仮想空間内に定義される複数の粒子を含む系を使用して流れをシミュレートする際、運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施すことを含む。
【００１１】
この態様によると、系の粒子を流れの擾乱源と関連付けられた粒子とそうでない粒子とに分けることができる。
【００１２】
本発明の別の態様は、解析装置である。この装置は、複数の粒子を含む系を仮想空間内に定義する粒子定義部と、粒子定義部によって定義された系に、流れを表すよう制約を課す制約付与部と、制約付与部によって制約が課された系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部と、を備える。数値演算部は流れの擾乱源と関連付けられた粒子を特定する。本解析装置はさらに、特定された粒子に、そうでない粒子とは異なる処理を施す処理部を備える。
【００１３】
なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。
【発明の効果】
【００１４】
本発明によれば、流れのシミュレーションにおいて、より容易に粒子の動きを把握することができる。
【図面の簡単な説明】
【００１５】
【図１】実施の形態に係る解析装置の機能および構成を示すブロック図である。
【図２】図１の粒子データ保持部の一例を示すデータ構造図である。
【図３】図１の解析装置における一連の処理の一例を示すフローチャートである。
【図４】図１の解析装置によるシミュレーション結果の一例を示す代表画面図である。
【図５】計算例に係る計算条件を示す模式図である。
【図６】計算例に係る計算条件の基で行われたシミュレーションの結果の時間変化を示す模式図である。
【発明を実施するための形態】
【００１６】
以下、各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。
【００１７】
実施の形態に係る解析装置は、仮想空間内に複数の粒子を含む系を定義し、その系を使用して現実世界の流れをシミュレートする。この際、解析装置は、運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施す。例えば、解析装置は、関連づけられた粒子をそうでない粒子とは異なる態様でディスプレイに表示させる。これにより、流体解析において、粒子の動きを把握しやすくなる。
【００１８】
図１は、実施の形態に係る解析装置１００の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。
【００１９】
解析装置１００は１、２または３次元の仮想空間内に流れの擾乱源を設定する。解析装置１００はそのように擾乱源が設定された仮想空間内に複数の粒子を配置し、それらの粒子に流れが表されるように制約を課す。解析装置１００は、粒子の運動方程式を数値的に演算することにより、擾乱源を含む仮想空間内での流れをシミュレートする。特に本実施の形態では、解析装置１００は、仮想空間内に渦の発生源を設定し、流れによる渦の発生をシミュレートする。
【００２０】
なお、本実施の形態ではＭＤ法に倣って粒子系を解析する場合について説明するが、繰り込み群分子動力学法（Renormalized Molecular Dynamics）やＤＥＭ（Distinct Element Method）やＳＰＨ（Smoothed Particle Hydrodynamics）やＭＰＳ（Moving Particle Semi-implicit）などの他の粒子法に倣って粒子系を解析する場合にも、本実施の形態に係る技術的思想を適用できることは本明細書に触れた当業者には明らかである。
【００２１】
解析装置１００は入力装置１０２およびディスプレイ１０４と接続される。入力装置１０２は、解析装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。
【００２２】
解析装置１００は、擾乱源定義部１０６と、粒子定義部１０８と、制約付与部１１０と、数値演算部１２０と、物理量演算部１３２と、表示制御部１１８と、擾乱源データ保持部１１２と、粒子データ保持部１１４と、表示データ保持部１１６と、を備える。
【００２３】
擾乱源定義部１０６は、入力装置１０２を介してユーザから取得する入力情報に基づき、仮想空間内に渦の発生源を定義する。この定義の仕方には少なくとも以下の２通りがある。
定義１．擾乱源定義部１０６は、渦の発生源を仮想空間内の領域（以下、発生源領域と称す）として定義する。以下の演算において、粒子が発生源領域に侵入するとその粒子は強い斥力を受ける。
定義２．擾乱源定義部１０６は、渦の発生源の表面上に、初期位置とバネでつながれた粒子（以下、発生源粒子と称す）を配置する。
いずれの定義が採用されるにせよ、擾乱源定義部１０６は、定義された渦の発生源のデータを擾乱源データ保持部１１２に登録する。
【００２４】
擾乱源データ保持部１１２は、渦の発生源のデータを保持する。渦の発生源の定義として定義１が採用される場合は、擾乱源データ保持部１１２は発生源領域を規定するパラメータや関数と、斥力についてのパラメータや関数と、を保持する。渦の発生源の定義として定義２が採用される場合は、擾乱源データ保持部１１２は発生源粒子を特定するＩＤと、その粒子の位置と、その粒子の速度と、を対応付けて保持する。
【００２５】
粒子定義部１０８は、入力装置１０２を介してユーザから取得する入力情報に基づき、仮想空間内にＮ（Ｎは自然数）個の粒子からなる粒子系Ｓを定義する。粒子定義部１０８は、仮想空間内の渦の発生源を除く領域にＮ個の粒子を配置する。粒子定義部１０８は、配置された粒子を特定する粒子ＩＤと、その粒子の位置と、を対応付けて粒子データ保持部１１４に登録する。
なお、粒子を現実世界の原子または分子に対応させてもよい。
【００２６】
制約付与部１１０は、入力装置１０２を介してユーザから取得する入力情報に基づき、粒子定義部１０８によって定義された粒子系Ｓに、流れを表すよう制約を課す。制約付与部１１０は、粒子系Ｓが全体として所定の向き（以下、流れの向きと称す）に移動するよう各粒子に初期速度を付与する。例えば制約付与部１１０は、粒子系Ｓの全ての粒子に、同じ向き、同じ大きさの初期速度を付与する。制約付与部１１０は、粒子ＩＤと、付与された初期速度と、を対応付けて粒子データ保持部１１４に登録する。
【００２７】
なお、制約付与部１１０は、粒子系Ｓの全ての粒子に、同じ向きであるが異なる大きさの初期速度を付与してもよい。あるいはまた、制約付与部１１０は、粒子の初期速度を粒子に亘って平均した平均初期速度（速度ベクトル）が実質的に流れの向きを向くように、各粒子に初期速度を付与してもよい。あるいはまた、制約付与部１１０は、初期速度（速度ベクトル）が渦の発生源を向くように、各粒子に初期速度を付与してもよい。
【００２８】
以下では粒子系Ｓの粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。
【００２９】
数値演算部１２０は、制約付与部１１０によって制約が課された粒子系Ｓの各粒子の運動を支配する支配方程式を数値的に演算する。特に数値演算部１２０は、離散化された粒子の運動方程式にしたがった繰り返し演算を行う。数値演算部１２０は繰り返し演算の過程で渦の発生源と関連付けられた粒子を特定する。
数値演算部１２０は、力演算部１２２と、粒子状態演算部１２４と、状態更新部１２６と、終了条件判定部１２８と、マーキング部１３０と、を含む。
【００３０】
力演算部１２２は粒子データ保持部１１４に保持される粒子系Ｓのデータおよび擾乱源データ保持部１１２に保持される渦の発生源のデータを参照し、粒子系Ｓの各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。力演算部１２２は、粒子系Ｓのｉ番目（１≦ｉ≦Ｎ）の粒子について、そのｉ番目の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。力演算部１２２は、各近接粒子について、その近接粒子とｉ番目の粒子との間のポテンシャルエネルギ関数およびその近接粒子とｉ番目の粒子との距離に基づいて、その近接粒子がｉ番目の粒子に及ぼす力を演算する。特に力演算部１２２は、その近接粒子とｉ番目の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。力演算部１２２は、近接粒子がｉ番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、ｉ番目の粒子に働く力を算出する。
【００３１】
渦の発生源の定義として定義１が採用される場合、力演算部１２２は、粒子系Ｓのｉ番目の粒子について、その粒子が発生源領域内にあるか否かを判定する。発生源領域内にあると判定された場合、力演算部１２２は、ｉ番目の粒子に働く力に渦の発生源起因の斥力を加算する。
渦の発生源の定義として定義２が採用される場合、力演算部１２２は、発生源粒子も含めてｉ番目の粒子の近接粒子を決定する。
【００３２】
粒子状態演算部１２４は粒子データ保持部１１４に保持される粒子系Ｓのデータを参照し、粒子系Ｓの各粒子について、離散化された粒子の運動方程式に力演算部１２２によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。
【００３３】
粒子状態演算部１２４は、力演算部１２２によって演算された力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、粒子系Ｓのｉ番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき所定の微小な時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって演算された力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。
【００３４】
粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、粒子系Ｓのｉ番目の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。
【００３５】
状態更新部１２６は、粒子データ保持部１１４に保持される粒子系Ｓの各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。状態更新部１２６はまた、更新後の粒子データ保持部１１４に保持される粒子のデータを、繰り返し演算のサイクル数に応じた情報と共に表示データ保持部１１６に登録する。繰り返し演算のサイクル数に応じた情報は、粒子系Ｓの時間発展における時刻に対応する。
表示データ保持部１１６は、粒子系Ｓの時間発展における時刻に対応する情報と、粒子系Ｓの粒子のデータと、を対応付けて保持する。
【００３６】
終了条件判定部１２８は、数値演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、外部から終了の指示を受け付けたことである。終了条件判定部１２８は、終了条件が満たされる場合、数値演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を力演算部１２２に戻す。すると力演算部１２２は、状態更新部１２６によって更新された粒子の位置で再び力を演算する。
【００３７】
マーキング部１３０は、以下の２つのモードのいずれかにより、渦の発生源と関連付けられた粒子をマーキングする。
第１モード．マーキング部１３０は、繰り返し演算の過程で渦の発生源の境界層に侵入した粒子をマーキングする。境界層については例えば非特許文献４や非特許文献５に詳しく説明されているので、ここでは以下に簡単に説明するに止め詳述しない。
【００３８】
非特許文献４、非特許文献５によると、擾乱源の代表長さｌと境界層の厚さδとの比はレイノルズ数Ｒｅの−０．５乗に比例する。
【数１】
例えば平板での境界層の厚さδは、
【数２】
である。
【００３９】
マーキング部１３０は、数値演算部１２０における繰り返し演算の１サイクル内の任意のタイミングで粒子のマーキングを行う。マーキング部１３０は、状態更新部１２６によって粒子データ保持部１１４が更新されると粒子のマーキングを行ってもよいし、力演算部１２２による演算の前に粒子のマーキングを行ってもよい。
【００４０】
マーキング部１３０は、擾乱源データ保持部１１２に保持される渦の発生源のデータに基づき境界層を決定する。マーキング部１３０は、粒子データ保持部１１４を参照し、粒子系Ｓの各粒子について、その粒子が境界層に入ったか否かを判定する。入ったと判定された場合、マーキング部１３０は、粒子データ保持部１１４においてその粒子の粒子ＩＤに対応付けられたマーキングフラグを立てる。マーキング部１３０は、マーキングフラグが既に立てられている場合は立てられたままとする。マーキングフラグは「０」または「１」の値をとり、「１」のときはマーキングフラグが立てられた状態であり、「０」のときはマーキングフラグが解除された状態である。
【００４１】
第２モード．マーキング部１３０は、繰り返し演算の過程で渦の発生源と相互作用した粒子をマーキングする。
渦の発生源の定義として定義１が採用される場合、マーキング部１３０は粒子データ保持部１１４および擾乱源データ保持部１１２を参照し、粒子系Ｓの各粒子について、その粒子が発生源と接触したかすなわちその粒子が発生源領域内に入ったか否かを判定する。入ったと判定された場合、マーキング部１３０は、粒子データ保持部１１４においてその粒子の粒子ＩＤに対応付けられたマーキングフラグを立てる。
【００４２】
渦の発生源の定義として定義２が採用される場合、マーキング部１３０は粒子データ保持部１１４および擾乱源データ保持部１１２を参照し、粒子系Ｓの各粒子について、発生源粒子がその粒子に力を及ぼすか否かを判定する。特にマーキング部１３０は、粒子系Ｓの各粒子について、その粒子との距離がカットオフ距離よりも小さくなる発生源粒子が存在するか否かを判定する。そのような発生源粒子が存在すると判定された場合、マーキング部１３０は、粒子データ保持部１１４においてその粒子の粒子ＩＤに対応付けられたマーキングフラグを立てる。
なお、マーキング部１３０の第２モードにおける処理は、力演算部１２２の処理の一部とされてもよい。
【００４３】
解析装置１００では、数値演算部１２０において特定された粒子に、そうでない粒子とは異なる処理が施される。例えば、特定された粒子とそうでない粒子とで別々に物理量が演算されてもよい。あるいはまた、特定された粒子とそうでない粒子とが異なる態様でディスプレイ１０４に表示されてもよい。
【００４４】
物理量演算部１３２は、粒子データ保持部１１４に保持される粒子系Ｓのデータに基づき粒子系Ｓの各種物理量、例えば温度や圧力や応力や渦度などを演算する。この際、物理量演算部１３２は、マーキングフラグが立っている粒子とそうでない粒子とで別々に演算を行ってもよい。あるいはまた、物理量演算部１３２は、マーキングフラグが立っている粒子についてのみ物理量を演算してもよく、また逆にそうでない粒子についてのみ物理量を演算してもよい。
【００４５】
表示制御部１１８は、表示データ保持部１１６に保持される粒子系Ｓの各粒子の位置、速度、マーキングフラグに基づき、ディスプレイ１０４に粒子系Ｓの時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。特に表示制御部１１８は、マーキングフラグが立っている粒子をそうでない粒子とは異なる態様でディスプレイ１０４に表示させる。
【００４６】
例えば、表示制御部１１８は、マーキングフラグが立っている粒子のみをディスプレイ１０４に表示させてもよいし、マーキングフラグが解除されている粒子のみをディスプレイ１０４に表示させてもよい。この場合、表示させるべき粒子の数を低減できるので、粒子の表示が表示用ハードウエアの性能に制限されにくくなる。あるいはまた、表示制御部１１８は、マーキングフラグが立っている粒子の色をそうでない粒子の色と異ならせてもよい。あるいはまた、表示制御部１１８は、物理量演算部１３２によって演算される渦度によるコンターを付加的に表示させてもよい。この場合、流れにより生じる渦の解析がより容易となる。
【００４７】
図２は、粒子データ保持部１１４の一例を示すデータ構造図である。粒子データ保持部１１４は、粒子ＩＤと、粒子の位置と、粒子の速度と、マーキングフラグと、を対応付けて保持する。
【００４８】
上述の実施の形態において、保持部の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。
【００４９】
以上の構成による解析装置１００の動作を説明する。
図３は、解析装置１００における一連の処理の一例を示すフローチャートである。図３はマーキング部１３０において第１モードが使用される場合を示す。擾乱源定義部１０６、粒子定義部１０８および制約付与部１１０は、渦の発生源のモデルと粒子系Ｓと各粒子の初期位置、初期速度とを含む初期条件を設定する（Ｓ２０２）。力演算部１２２は、粒子間の距離から粒子に働く力を演算する（Ｓ２０４）。粒子状態演算部１２４は、演算された力を含む粒子の運動方程式から速度を演算する（Ｓ２０６）。粒子状態演算部１２４は、演算された速度から粒子の位置を算出する（Ｓ２０８）。状態更新部１２６は、粒子データ保持部１１４に保持される粒子の位置、速度を算出された位置、速度で更新する（Ｓ２１０）。マーキング部１３０は、更新後の粒子データ保持部１１４を参照し、新たに境界層に侵入した粒子を特定する（Ｓ２１２）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２１４）。終了条件が満たされない場合（Ｓ２１４のＮ）、処理はＳ２０４に戻される。終了条件が満たされる場合（Ｓ２１４のＹ）、ディスプレイ１０４は特定された粒子を識別可能に表示する（Ｓ２１６）。
【００５０】
本実施の形態に係る解析装置１００によると、繰り返し演算の過程で渦の発生源と関連付けられた粒子が特定され、そのように特定された粒子はそうでない粒子とは異なる態様で処理される。したがって、渦の発生源が流れに及ぼす作用をより正確にシミュレートできる。
【００５１】
また、渦の発生源と関連付けられた粒子はそうでない粒子とは異なる態様でディスプレイ１０４に表示されるので、渦をより正確に可視化することができる。また、渦の位置や大きさと同時に流体粒子の動きもより容易に把握することができる。
【００５２】
流れを可視化する実験の手法としては、懸濁法や電解沈殿法（例えば、非特許文献６参照）などが知られている。本発明者は、流れのシミュレーションにおいても電解沈殿法のような可視化ができないかという着想に基づき、本実施の形態に係る解析装置１００を創作した。電解沈殿法では、例えば水を電気分解するとき陽極近傍に生成される白色の沈殿物をトレーサとして利用する。そこで本発明者は、流れのシミュレーションにおいて、渦の発生源の境界層、または、その発生源との相互作用の有無、をマーキングの基準とすることにより、電解沈殿法に倣って渦をより効果的に可視化できることに想到した。
【００５３】
なお、懸濁法に倣い、粒子系Ｓに粒子系Ｓの粒子とは性質が異なる粒子を最初から加えておき、その異粒子を追跡することで流れをシミュレートすることも考えられる。しかしながら、本実施の形態に係る手法は少なくとも、運動の過程で渦の発生源と関連付けられた粒子を特定する点で、上記の懸濁法に倣った手法とは異なる。本実施の形態に係る手法では、粒子系Ｓに異粒子を加えることは行われない。したがって本実施の形態に係る手法では、異粒子の付加により粒子系Ｓが乱されることはないので、流れをより正確にシミュレートできる。また、異粒子を加える手法では渦とは関係のない部分でも異粒子が観察されるのに対し、本実施の形態に係る手法では渦の発生源と関連付けられた粒子が特定される。したがって、可視化したときには前者よりも後者のほうがより明確に渦の状態を示す。
【００５４】
図４は、解析装置１００によるシミュレーション結果の一例を示す代表画面図である。図４ではマーキング部１３０によってマーキングされた粒子のみが色つきで示されている。この図４と、非特許文献６に示される電解沈殿法により得られた円柱後方のカルマン渦列の図とを比較すると、本実施の形態に係る手法が渦をより良くシミュレートしていることが分かる。
【００５５】
本発明者は、本実施の形態に係る手法の効果を確認するため、ＭＤ法によりシリンダ周りの流れをシミュレートし、本実施の形態に係る可視化方法を適用した。
図５は、計算例に係る計算条件を示す模式図である。２次元の仮想空間が設定された。仮想空間はｘ方向３４６４×１０−１０（ｍ）、ｙ方向３４６４×１０−１０（ｍ）とし、粒子系Ｓの粒子数を約８０万個、粒子の質量を６．６３×１０−２６（ｋｇ）とした。シリンダの直径は５０５×１０−１０（ｍ）とした。流入の温度を８５（Ｋ）、初期速度をｘ軸の正の向きに３００（ｍ／ｓ）とし、さらに粒子に加速度７．０８×１０１１（ｍ／ｓ２）を加えた。この加速度はシリンダを配置したことでの圧力損失によって粒子がｘ軸の負の向きに戻るのを防ぐために付加されている。粒子間の相互作用は以下に示す式（１）のようにＭＤ法でよく用いられるレナードジョーンズポテンシャルの斥力部分のみを用いた。ポテンシャルパラメータは以下に示す式（２）のように設定した。
【数３】
さらに渦の発生源の定義として定義２を採用し、シリンダの表面上(本計算例では２次元のため円周上)に初期位置とバネでつながれた発生源粒子を配置した。そしてシリンダ表面に配置した発生源粒子にはランジュバン法による温度制御（非特許文献７参照）を施した。シリンダの温度は８５（Ｋ）に設定した。
【００５６】
粒子系Ｓの粒子の粘度は３．６４×１０−４（Ｐａ・ｓ）（非特許文献８参照）、密度は１．１８×１０３（ｋｇ／ｍ３）である。代表長さはシリンダの直径であり、５０５×１０−１０（ｍ）である。さらに流速はシリンダ前方の平均速度より４２６（ｍ／ｓ）である。よってこの系のレイノルズ数は７０．１である。
【００５７】
図６は、計算例に係る計算条件の基で行われたシミュレーションの結果の時間変化を示す模式図である。図６では図４と同様に、マーキングされた粒子のみが色つきで示されている。図６に示される通り、本実施の形態に係る手法によって、ＭＤ法において電解沈殿法のような可視化が可能であることが確認された。さらに、従来の手法による可視化と比較すると、本実施の形態に係る手法ではより容易に、渦の位置や大きさと同時に粒子の動きを把握することができる。すなわち、本実施の形態に係る手法は流体解析の可視化に有用な手法である。
【００５８】
以上、実施の形態に係る解析装置１００の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。
【００５９】
実施の形態では、数値演算部１２０において粒子の位置と速度の両方を演算する場合について説明したが、これに限られない。例えば、数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法もあり、本実施の形態に係る技術的思想をそのような手法に適用してもよい。
【００６０】
実施の形態では、マーキング部１３０は第１モードにおいて、繰り返し演算の過程で渦の発生源の境界層に侵入した粒子をマーキングする場合について説明したが、これに限られない。例えば境界層よりも任意に薄い層を定義し、マーキング部は繰り返し演算の過程でそのように定義された薄い層に侵入した粒子をマーキングしてもよい。
【符号の説明】
【００６１】
１００ 解析装置、 １０２ 入力装置、 １０４ ディスプレイ、 １０６ 擾乱源定義部、 １０８ 粒子定義部、 １１０ 制約付与部、 １１４ 粒子データ保持部、 １１６ 表示データ保持部、 １１８ 表示制御部、 １２０ 数値演算部、 １３２ 物理量演算部。
【技術分野】
【０００１】
本発明は、粒子系を使用したシミュレーション方法および粒子系を解析する解析装置に関する。
【背景技術】
【０００２】
従来、古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）に基づくシミュレーションが知られている。ＭＤ法は物性をポテンシャルエネルギ関数で与えるのでモデル化の少ない手法である。
【０００３】
非特許文献１、２は、ＭＤ法による流体解析について報告している。非特許文献１、２では、ＭＤ法を用いてシリンダ周りの流れの計算が行われ、速度ベクトルの分布によってカルマン渦の可視化が試みられている。より具体的には、領域が分割され、分割された領域ごとの平均速度の向きが矢印で示されている。
【０００４】
また、ＭＤ法以外の粒子法（例えば、特許文献１、非特許文献３参照）では、速度分布以外にも、圧力分布や密度分布などによって渦を可視化することが試みられている。
有限要素法でも、速度分布や圧力分布以外にも渦度の分布から渦の確認が試みられている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００６−２８５８６６号公報
【非特許文献】
【０００６】
【非特許文献１】Ｄ．Ｃ．Ｒａｐａｐｏｒｔ，Ｅ．Ｃｌｅｍｅｎｔｉ、Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．、１９８６、ｖｏｌ５７、ｎｏ６、ｐ．６９５−６９８
【非特許文献２】Ｄ．Ｃ．Ｒａｐａｐｏｒｔ、Ｐｈｙｓ．Ｒｅｖ．Ａ、１９８７、ｖｏｌ３６、ｎｏ７、ｐ．３２８８−３２９９
【非特許文献３】越塚誠一、粒子法（計算力学レクチャーシリーズ５）、日本計算工学会編、丸善、２００５年
【非特許文献４】池畑光尚、船舶海洋工学のための流体力学入門、船舶技術協会、１９９３年、ｐ．１３０
【非特許文献５】白倉昌明、大橋秀雄、標準機械工学講座２３ 流体力学（２）、コロナ社、１９６９年、第２章
【非特許文献６】日本機械学会編、流れ−写真集、丸善、１９８４年、ｐ．４
【非特許文献７】丸山茂夫、井上和洋、「第１３回計算力学講演会講演論文集」、２０００年、ｐ．３８９−３９０
【非特許文献８】Ｄ．Ｌｅｖｅｓｑｕｅ，Ｌ．Ｖｅｒｌｅｔ，Ｊ．Ｋｕｒｋｉｊａｒｖｉ、Ｐｈｙｓ．Ｒｅｖ．Ａ、１９７３、ｖｏｌ７、ｎｏ５、ｐ．１６９０−１７００
【発明の概要】
【発明が解決しようとする課題】
【０００７】
上記の従来手法はいずれも、仮想空間内の単位領域ごとに速度や圧力などの物理量を計算し、計算された物理量をグラデーション等により見やすくしてユーザに提供することにより渦の可視化を目指している。このような手法は定常状態の解析には向いているかもしれないが、渦のように生来的に動的な現象をシミュレートする場合、流体粒子の動きを把握することが困難であり不向きである。
【０００８】
このような課題は渦のシミュレーションに限らず一般的な流れの擾乱のシミュレーションでも起こりうる。
【０００９】
本発明はこうした課題に鑑みてなされたものであり、その目的は、流れのシミュレーションにおいて、より容易に粒子の動きを把握することができる解析技術の提供にある。
【課題を解決するための手段】
【００１０】
本発明のある態様は、シミュレーション方法に関する。このシミュレーション方法は、仮想空間内に定義される複数の粒子を含む系を使用して流れをシミュレートする際、運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施すことを含む。
【００１１】
この態様によると、系の粒子を流れの擾乱源と関連付けられた粒子とそうでない粒子とに分けることができる。
【００１２】
本発明の別の態様は、解析装置である。この装置は、複数の粒子を含む系を仮想空間内に定義する粒子定義部と、粒子定義部によって定義された系に、流れを表すよう制約を課す制約付与部と、制約付与部によって制約が課された系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部と、を備える。数値演算部は流れの擾乱源と関連付けられた粒子を特定する。本解析装置はさらに、特定された粒子に、そうでない粒子とは異なる処理を施す処理部を備える。
【００１３】
なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。
【発明の効果】
【００１４】
本発明によれば、流れのシミュレーションにおいて、より容易に粒子の動きを把握することができる。
【図面の簡単な説明】
【００１５】
【図１】実施の形態に係る解析装置の機能および構成を示すブロック図である。
【図２】図１の粒子データ保持部の一例を示すデータ構造図である。
【図３】図１の解析装置における一連の処理の一例を示すフローチャートである。
【図４】図１の解析装置によるシミュレーション結果の一例を示す代表画面図である。
【図５】計算例に係る計算条件を示す模式図である。
【図６】計算例に係る計算条件の基で行われたシミュレーションの結果の時間変化を示す模式図である。
【発明を実施するための形態】
【００１６】
以下、各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。
【００１７】
実施の形態に係る解析装置は、仮想空間内に複数の粒子を含む系を定義し、その系を使用して現実世界の流れをシミュレートする。この際、解析装置は、運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施す。例えば、解析装置は、関連づけられた粒子をそうでない粒子とは異なる態様でディスプレイに表示させる。これにより、流体解析において、粒子の動きを把握しやすくなる。
【００１８】
図１は、実施の形態に係る解析装置１００の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。
【００１９】
解析装置１００は１、２または３次元の仮想空間内に流れの擾乱源を設定する。解析装置１００はそのように擾乱源が設定された仮想空間内に複数の粒子を配置し、それらの粒子に流れが表されるように制約を課す。解析装置１００は、粒子の運動方程式を数値的に演算することにより、擾乱源を含む仮想空間内での流れをシミュレートする。特に本実施の形態では、解析装置１００は、仮想空間内に渦の発生源を設定し、流れによる渦の発生をシミュレートする。
【００２０】
なお、本実施の形態ではＭＤ法に倣って粒子系を解析する場合について説明するが、繰り込み群分子動力学法（Renormalized Molecular Dynamics）やＤＥＭ（Distinct Element Method）やＳＰＨ（Smoothed Particle Hydrodynamics）やＭＰＳ（Moving Particle Semi-implicit）などの他の粒子法に倣って粒子系を解析する場合にも、本実施の形態に係る技術的思想を適用できることは本明細書に触れた当業者には明らかである。
【００２１】
解析装置１００は入力装置１０２およびディスプレイ１０４と接続される。入力装置１０２は、解析装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。
【００２２】
解析装置１００は、擾乱源定義部１０６と、粒子定義部１０８と、制約付与部１１０と、数値演算部１２０と、物理量演算部１３２と、表示制御部１１８と、擾乱源データ保持部１１２と、粒子データ保持部１１４と、表示データ保持部１１６と、を備える。
【００２３】
擾乱源定義部１０６は、入力装置１０２を介してユーザから取得する入力情報に基づき、仮想空間内に渦の発生源を定義する。この定義の仕方には少なくとも以下の２通りがある。
定義１．擾乱源定義部１０６は、渦の発生源を仮想空間内の領域（以下、発生源領域と称す）として定義する。以下の演算において、粒子が発生源領域に侵入するとその粒子は強い斥力を受ける。
定義２．擾乱源定義部１０６は、渦の発生源の表面上に、初期位置とバネでつながれた粒子（以下、発生源粒子と称す）を配置する。
いずれの定義が採用されるにせよ、擾乱源定義部１０６は、定義された渦の発生源のデータを擾乱源データ保持部１１２に登録する。
【００２４】
擾乱源データ保持部１１２は、渦の発生源のデータを保持する。渦の発生源の定義として定義１が採用される場合は、擾乱源データ保持部１１２は発生源領域を規定するパラメータや関数と、斥力についてのパラメータや関数と、を保持する。渦の発生源の定義として定義２が採用される場合は、擾乱源データ保持部１１２は発生源粒子を特定するＩＤと、その粒子の位置と、その粒子の速度と、を対応付けて保持する。
【００２５】
粒子定義部１０８は、入力装置１０２を介してユーザから取得する入力情報に基づき、仮想空間内にＮ（Ｎは自然数）個の粒子からなる粒子系Ｓを定義する。粒子定義部１０８は、仮想空間内の渦の発生源を除く領域にＮ個の粒子を配置する。粒子定義部１０８は、配置された粒子を特定する粒子ＩＤと、その粒子の位置と、を対応付けて粒子データ保持部１１４に登録する。
なお、粒子を現実世界の原子または分子に対応させてもよい。
【００２６】
制約付与部１１０は、入力装置１０２を介してユーザから取得する入力情報に基づき、粒子定義部１０８によって定義された粒子系Ｓに、流れを表すよう制約を課す。制約付与部１１０は、粒子系Ｓが全体として所定の向き（以下、流れの向きと称す）に移動するよう各粒子に初期速度を付与する。例えば制約付与部１１０は、粒子系Ｓの全ての粒子に、同じ向き、同じ大きさの初期速度を付与する。制約付与部１１０は、粒子ＩＤと、付与された初期速度と、を対応付けて粒子データ保持部１１４に登録する。
【００２７】
なお、制約付与部１１０は、粒子系Ｓの全ての粒子に、同じ向きであるが異なる大きさの初期速度を付与してもよい。あるいはまた、制約付与部１１０は、粒子の初期速度を粒子に亘って平均した平均初期速度（速度ベクトル）が実質的に流れの向きを向くように、各粒子に初期速度を付与してもよい。あるいはまた、制約付与部１１０は、初期速度（速度ベクトル）が渦の発生源を向くように、各粒子に初期速度を付与してもよい。
【００２８】
以下では粒子系Ｓの粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。
【００２９】
数値演算部１２０は、制約付与部１１０によって制約が課された粒子系Ｓの各粒子の運動を支配する支配方程式を数値的に演算する。特に数値演算部１２０は、離散化された粒子の運動方程式にしたがった繰り返し演算を行う。数値演算部１２０は繰り返し演算の過程で渦の発生源と関連付けられた粒子を特定する。
数値演算部１２０は、力演算部１２２と、粒子状態演算部１２４と、状態更新部１２６と、終了条件判定部１２８と、マーキング部１３０と、を含む。
【００３０】
力演算部１２２は粒子データ保持部１１４に保持される粒子系Ｓのデータおよび擾乱源データ保持部１１２に保持される渦の発生源のデータを参照し、粒子系Ｓの各粒子について、粒子間の距離に基づきその粒子に働く力を演算する。力演算部１２２は、粒子系Ｓのｉ番目（１≦ｉ≦Ｎ）の粒子について、そのｉ番目の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。力演算部１２２は、各近接粒子について、その近接粒子とｉ番目の粒子との間のポテンシャルエネルギ関数およびその近接粒子とｉ番目の粒子との距離に基づいて、その近接粒子がｉ番目の粒子に及ぼす力を演算する。特に力演算部１２２は、その近接粒子とｉ番目の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。力演算部１２２は、近接粒子がｉ番目の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、ｉ番目の粒子に働く力を算出する。
【００３１】
渦の発生源の定義として定義１が採用される場合、力演算部１２２は、粒子系Ｓのｉ番目の粒子について、その粒子が発生源領域内にあるか否かを判定する。発生源領域内にあると判定された場合、力演算部１２２は、ｉ番目の粒子に働く力に渦の発生源起因の斥力を加算する。
渦の発生源の定義として定義２が採用される場合、力演算部１２２は、発生源粒子も含めてｉ番目の粒子の近接粒子を決定する。
【００３２】
粒子状態演算部１２４は粒子データ保持部１１４に保持される粒子系Ｓのデータを参照し、粒子系Ｓの各粒子について、離散化された粒子の運動方程式に力演算部１２２によって演算された力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。
【００３３】
粒子状態演算部１２４は、力演算部１２２によって演算された力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、粒子系Ｓのｉ番目の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき所定の微小な時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって演算された力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。
【００３４】
粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、粒子系Ｓのｉ番目の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。
【００３５】
状態更新部１２６は、粒子データ保持部１１４に保持される粒子系Ｓの各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。状態更新部１２６はまた、更新後の粒子データ保持部１１４に保持される粒子のデータを、繰り返し演算のサイクル数に応じた情報と共に表示データ保持部１１６に登録する。繰り返し演算のサイクル数に応じた情報は、粒子系Ｓの時間発展における時刻に対応する。
表示データ保持部１１６は、粒子系Ｓの時間発展における時刻に対応する情報と、粒子系Ｓの粒子のデータと、を対応付けて保持する。
【００３６】
終了条件判定部１２８は、数値演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、外部から終了の指示を受け付けたことである。終了条件判定部１２８は、終了条件が満たされる場合、数値演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を力演算部１２２に戻す。すると力演算部１２２は、状態更新部１２６によって更新された粒子の位置で再び力を演算する。
【００３７】
マーキング部１３０は、以下の２つのモードのいずれかにより、渦の発生源と関連付けられた粒子をマーキングする。
第１モード．マーキング部１３０は、繰り返し演算の過程で渦の発生源の境界層に侵入した粒子をマーキングする。境界層については例えば非特許文献４や非特許文献５に詳しく説明されているので、ここでは以下に簡単に説明するに止め詳述しない。
【００３８】
非特許文献４、非特許文献５によると、擾乱源の代表長さｌと境界層の厚さδとの比はレイノルズ数Ｒｅの−０．５乗に比例する。
【数１】
例えば平板での境界層の厚さδは、
【数２】
である。
【００３９】
マーキング部１３０は、数値演算部１２０における繰り返し演算の１サイクル内の任意のタイミングで粒子のマーキングを行う。マーキング部１３０は、状態更新部１２６によって粒子データ保持部１１４が更新されると粒子のマーキングを行ってもよいし、力演算部１２２による演算の前に粒子のマーキングを行ってもよい。
【００４０】
マーキング部１３０は、擾乱源データ保持部１１２に保持される渦の発生源のデータに基づき境界層を決定する。マーキング部１３０は、粒子データ保持部１１４を参照し、粒子系Ｓの各粒子について、その粒子が境界層に入ったか否かを判定する。入ったと判定された場合、マーキング部１３０は、粒子データ保持部１１４においてその粒子の粒子ＩＤに対応付けられたマーキングフラグを立てる。マーキング部１３０は、マーキングフラグが既に立てられている場合は立てられたままとする。マーキングフラグは「０」または「１」の値をとり、「１」のときはマーキングフラグが立てられた状態であり、「０」のときはマーキングフラグが解除された状態である。
【００４１】
第２モード．マーキング部１３０は、繰り返し演算の過程で渦の発生源と相互作用した粒子をマーキングする。
渦の発生源の定義として定義１が採用される場合、マーキング部１３０は粒子データ保持部１１４および擾乱源データ保持部１１２を参照し、粒子系Ｓの各粒子について、その粒子が発生源と接触したかすなわちその粒子が発生源領域内に入ったか否かを判定する。入ったと判定された場合、マーキング部１３０は、粒子データ保持部１１４においてその粒子の粒子ＩＤに対応付けられたマーキングフラグを立てる。
【００４２】
渦の発生源の定義として定義２が採用される場合、マーキング部１３０は粒子データ保持部１１４および擾乱源データ保持部１１２を参照し、粒子系Ｓの各粒子について、発生源粒子がその粒子に力を及ぼすか否かを判定する。特にマーキング部１３０は、粒子系Ｓの各粒子について、その粒子との距離がカットオフ距離よりも小さくなる発生源粒子が存在するか否かを判定する。そのような発生源粒子が存在すると判定された場合、マーキング部１３０は、粒子データ保持部１１４においてその粒子の粒子ＩＤに対応付けられたマーキングフラグを立てる。
なお、マーキング部１３０の第２モードにおける処理は、力演算部１２２の処理の一部とされてもよい。
【００４３】
解析装置１００では、数値演算部１２０において特定された粒子に、そうでない粒子とは異なる処理が施される。例えば、特定された粒子とそうでない粒子とで別々に物理量が演算されてもよい。あるいはまた、特定された粒子とそうでない粒子とが異なる態様でディスプレイ１０４に表示されてもよい。
【００４４】
物理量演算部１３２は、粒子データ保持部１１４に保持される粒子系Ｓのデータに基づき粒子系Ｓの各種物理量、例えば温度や圧力や応力や渦度などを演算する。この際、物理量演算部１３２は、マーキングフラグが立っている粒子とそうでない粒子とで別々に演算を行ってもよい。あるいはまた、物理量演算部１３２は、マーキングフラグが立っている粒子についてのみ物理量を演算してもよく、また逆にそうでない粒子についてのみ物理量を演算してもよい。
【００４５】
表示制御部１１８は、表示データ保持部１１６に保持される粒子系Ｓの各粒子の位置、速度、マーキングフラグに基づき、ディスプレイ１０４に粒子系Ｓの時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。特に表示制御部１１８は、マーキングフラグが立っている粒子をそうでない粒子とは異なる態様でディスプレイ１０４に表示させる。
【００４６】
例えば、表示制御部１１８は、マーキングフラグが立っている粒子のみをディスプレイ１０４に表示させてもよいし、マーキングフラグが解除されている粒子のみをディスプレイ１０４に表示させてもよい。この場合、表示させるべき粒子の数を低減できるので、粒子の表示が表示用ハードウエアの性能に制限されにくくなる。あるいはまた、表示制御部１１８は、マーキングフラグが立っている粒子の色をそうでない粒子の色と異ならせてもよい。あるいはまた、表示制御部１１８は、物理量演算部１３２によって演算される渦度によるコンターを付加的に表示させてもよい。この場合、流れにより生じる渦の解析がより容易となる。
【００４７】
図２は、粒子データ保持部１１４の一例を示すデータ構造図である。粒子データ保持部１１４は、粒子ＩＤと、粒子の位置と、粒子の速度と、マーキングフラグと、を対応付けて保持する。
【００４８】
上述の実施の形態において、保持部の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。
【００４９】
以上の構成による解析装置１００の動作を説明する。
図３は、解析装置１００における一連の処理の一例を示すフローチャートである。図３はマーキング部１３０において第１モードが使用される場合を示す。擾乱源定義部１０６、粒子定義部１０８および制約付与部１１０は、渦の発生源のモデルと粒子系Ｓと各粒子の初期位置、初期速度とを含む初期条件を設定する（Ｓ２０２）。力演算部１２２は、粒子間の距離から粒子に働く力を演算する（Ｓ２０４）。粒子状態演算部１２４は、演算された力を含む粒子の運動方程式から速度を演算する（Ｓ２０６）。粒子状態演算部１２４は、演算された速度から粒子の位置を算出する（Ｓ２０８）。状態更新部１２６は、粒子データ保持部１１４に保持される粒子の位置、速度を算出された位置、速度で更新する（Ｓ２１０）。マーキング部１３０は、更新後の粒子データ保持部１１４を参照し、新たに境界層に侵入した粒子を特定する（Ｓ２１２）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２１４）。終了条件が満たされない場合（Ｓ２１４のＮ）、処理はＳ２０４に戻される。終了条件が満たされる場合（Ｓ２１４のＹ）、ディスプレイ１０４は特定された粒子を識別可能に表示する（Ｓ２１６）。
【００５０】
本実施の形態に係る解析装置１００によると、繰り返し演算の過程で渦の発生源と関連付けられた粒子が特定され、そのように特定された粒子はそうでない粒子とは異なる態様で処理される。したがって、渦の発生源が流れに及ぼす作用をより正確にシミュレートできる。
【００５１】
また、渦の発生源と関連付けられた粒子はそうでない粒子とは異なる態様でディスプレイ１０４に表示されるので、渦をより正確に可視化することができる。また、渦の位置や大きさと同時に流体粒子の動きもより容易に把握することができる。
【００５２】
流れを可視化する実験の手法としては、懸濁法や電解沈殿法（例えば、非特許文献６参照）などが知られている。本発明者は、流れのシミュレーションにおいても電解沈殿法のような可視化ができないかという着想に基づき、本実施の形態に係る解析装置１００を創作した。電解沈殿法では、例えば水を電気分解するとき陽極近傍に生成される白色の沈殿物をトレーサとして利用する。そこで本発明者は、流れのシミュレーションにおいて、渦の発生源の境界層、または、その発生源との相互作用の有無、をマーキングの基準とすることにより、電解沈殿法に倣って渦をより効果的に可視化できることに想到した。
【００５３】
なお、懸濁法に倣い、粒子系Ｓに粒子系Ｓの粒子とは性質が異なる粒子を最初から加えておき、その異粒子を追跡することで流れをシミュレートすることも考えられる。しかしながら、本実施の形態に係る手法は少なくとも、運動の過程で渦の発生源と関連付けられた粒子を特定する点で、上記の懸濁法に倣った手法とは異なる。本実施の形態に係る手法では、粒子系Ｓに異粒子を加えることは行われない。したがって本実施の形態に係る手法では、異粒子の付加により粒子系Ｓが乱されることはないので、流れをより正確にシミュレートできる。また、異粒子を加える手法では渦とは関係のない部分でも異粒子が観察されるのに対し、本実施の形態に係る手法では渦の発生源と関連付けられた粒子が特定される。したがって、可視化したときには前者よりも後者のほうがより明確に渦の状態を示す。
【００５４】
図４は、解析装置１００によるシミュレーション結果の一例を示す代表画面図である。図４ではマーキング部１３０によってマーキングされた粒子のみが色つきで示されている。この図４と、非特許文献６に示される電解沈殿法により得られた円柱後方のカルマン渦列の図とを比較すると、本実施の形態に係る手法が渦をより良くシミュレートしていることが分かる。
【００５５】
本発明者は、本実施の形態に係る手法の効果を確認するため、ＭＤ法によりシリンダ周りの流れをシミュレートし、本実施の形態に係る可視化方法を適用した。
図５は、計算例に係る計算条件を示す模式図である。２次元の仮想空間が設定された。仮想空間はｘ方向３４６４×１０−１０（ｍ）、ｙ方向３４６４×１０−１０（ｍ）とし、粒子系Ｓの粒子数を約８０万個、粒子の質量を６．６３×１０−２６（ｋｇ）とした。シリンダの直径は５０５×１０−１０（ｍ）とした。流入の温度を８５（Ｋ）、初期速度をｘ軸の正の向きに３００（ｍ／ｓ）とし、さらに粒子に加速度７．０８×１０１１（ｍ／ｓ２）を加えた。この加速度はシリンダを配置したことでの圧力損失によって粒子がｘ軸の負の向きに戻るのを防ぐために付加されている。粒子間の相互作用は以下に示す式（１）のようにＭＤ法でよく用いられるレナードジョーンズポテンシャルの斥力部分のみを用いた。ポテンシャルパラメータは以下に示す式（２）のように設定した。
【数３】
さらに渦の発生源の定義として定義２を採用し、シリンダの表面上(本計算例では２次元のため円周上)に初期位置とバネでつながれた発生源粒子を配置した。そしてシリンダ表面に配置した発生源粒子にはランジュバン法による温度制御（非特許文献７参照）を施した。シリンダの温度は８５（Ｋ）に設定した。
【００５６】
粒子系Ｓの粒子の粘度は３．６４×１０−４（Ｐａ・ｓ）（非特許文献８参照）、密度は１．１８×１０３（ｋｇ／ｍ３）である。代表長さはシリンダの直径であり、５０５×１０−１０（ｍ）である。さらに流速はシリンダ前方の平均速度より４２６（ｍ／ｓ）である。よってこの系のレイノルズ数は７０．１である。
【００５７】
図６は、計算例に係る計算条件の基で行われたシミュレーションの結果の時間変化を示す模式図である。図６では図４と同様に、マーキングされた粒子のみが色つきで示されている。図６に示される通り、本実施の形態に係る手法によって、ＭＤ法において電解沈殿法のような可視化が可能であることが確認された。さらに、従来の手法による可視化と比較すると、本実施の形態に係る手法ではより容易に、渦の位置や大きさと同時に粒子の動きを把握することができる。すなわち、本実施の形態に係る手法は流体解析の可視化に有用な手法である。
【００５８】
以上、実施の形態に係る解析装置１００の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。
【００５９】
実施の形態では、数値演算部１２０において粒子の位置と速度の両方を演算する場合について説明したが、これに限られない。例えば、数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法もあり、本実施の形態に係る技術的思想をそのような手法に適用してもよい。
【００６０】
実施の形態では、マーキング部１３０は第１モードにおいて、繰り返し演算の過程で渦の発生源の境界層に侵入した粒子をマーキングする場合について説明したが、これに限られない。例えば境界層よりも任意に薄い層を定義し、マーキング部は繰り返し演算の過程でそのように定義された薄い層に侵入した粒子をマーキングしてもよい。
【符号の説明】
【００６１】
１００ 解析装置、 １０２ 入力装置、 １０４ ディスプレイ、 １０６ 擾乱源定義部、 １０８ 粒子定義部、 １１０ 制約付与部、 １１４ 粒子データ保持部、 １１６ 表示データ保持部、 １１８ 表示制御部、 １２０ 数値演算部、 １３２ 物理量演算部。
【特許請求の範囲】
【請求項１】
仮想空間内に定義される複数の粒子を含む系を使用して流れをシミュレートする際、運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施すことを特徴とするシミュレーション方法。
【請求項２】
前記擾乱源の境界層に侵入した粒子に、そうでない粒子とは異なる処理を施すことを特徴とする請求項１に記載のシミュレーション方法。
【請求項３】
前記擾乱源と相互作用した粒子に、そうでない粒子とは異なる処理を施すことを特徴とする請求項１または２に記載のシミュレーション方法。
【請求項４】
所定の表示装置に、前記擾乱源と関連付けられた粒子をそうでない粒子とは異なる態様で表示させることを特徴とする請求項１から３のいずれかに記載のシミュレーション方法。
【請求項５】
複数の粒子を含む系を仮想空間内に定義する粒子定義部と、
前記粒子定義部によって定義された系に、流れを表すよう制約を課す制約付与部と、
前記制約付与部によって制約が課された系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部と、を備え、
前記数値演算部は流れの擾乱源と関連付けられた粒子を特定し、
本解析装置はさらに、
特定された粒子に、そうでない粒子とは異なる処理を施す処理部を備えることを特徴とする解析装置。
【請求項６】
仮想空間内に定義される複数の粒子を含む系を使用して流れをシミュレートする機能をコンピュータに実現させるコンピュータプログラムであって、
運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施す機能を前記コンピュータに実現させることを特徴とするコンピュータプログラム。
【請求項１】
仮想空間内に定義される複数の粒子を含む系を使用して流れをシミュレートする際、運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施すことを特徴とするシミュレーション方法。
【請求項２】
前記擾乱源の境界層に侵入した粒子に、そうでない粒子とは異なる処理を施すことを特徴とする請求項１に記載のシミュレーション方法。
【請求項３】
前記擾乱源と相互作用した粒子に、そうでない粒子とは異なる処理を施すことを特徴とする請求項１または２に記載のシミュレーション方法。
【請求項４】
所定の表示装置に、前記擾乱源と関連付けられた粒子をそうでない粒子とは異なる態様で表示させることを特徴とする請求項１から３のいずれかに記載のシミュレーション方法。
【請求項５】
複数の粒子を含む系を仮想空間内に定義する粒子定義部と、
前記粒子定義部によって定義された系に、流れを表すよう制約を課す制約付与部と、
前記制約付与部によって制約が課された系の各粒子の運動を支配する支配方程式を数値的に演算する数値演算部と、を備え、
前記数値演算部は流れの擾乱源と関連付けられた粒子を特定し、
本解析装置はさらに、
特定された粒子に、そうでない粒子とは異なる処理を施す処理部を備えることを特徴とする解析装置。
【請求項６】
仮想空間内に定義される複数の粒子を含む系を使用して流れをシミュレートする機能をコンピュータに実現させるコンピュータプログラムであって、
運動の過程で流れの擾乱源と関連付けられた粒子に、そうでない粒子とは異なる処理を施す機能を前記コンピュータに実現させることを特徴とするコンピュータプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図２】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２０１３−７３５４０（Ｐ２０１３−７３５４０Ａ）
【公開日】平成２５年４月２２日（２０１３．４．２２）
【国際特許分類】
【出願番号】特願２０１１−２１３７１２（Ｐ２０１１−２１３７１２）
【出願日】平成２３年９月２９日（２０１１．９．２９）
【出願人】（０００００２１０７）住友重機械工業株式会社 (2,241)
【公開日】平成２５年４月２２日（２０１３．４．２２）
【国際特許分類】
【出願日】平成２３年９月２９日（２０１１．９．２９）
【出願人】（０００００２１０７）住友重機械工業株式会社 (2,241)
[ Back to top ]