JP6458501B2 - シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置 - Google Patents

シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置 Download PDF

Info

Publication number
JP6458501B2
JP6458501B2 JP2015000797A JP2015000797A JP6458501B2 JP 6458501 B2 JP6458501 B2 JP 6458501B2 JP 2015000797 A JP2015000797 A JP 2015000797A JP 2015000797 A JP2015000797 A JP 2015000797A JP 6458501 B2 JP6458501 B2 JP 6458501B2
Authority
JP
Japan
Prior art keywords
particles
region
particle
boundary
repulsive force
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2015000797A
Other languages
English (en)
Other versions
JP2016125934A (ja
Inventor
多聞 諏訪
多聞 諏訪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujitsu Ltd
Original Assignee
Fujitsu Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujitsu Ltd filed Critical Fujitsu Ltd
Priority to JP2015000797A priority Critical patent/JP6458501B2/ja
Priority to US14/950,373 priority patent/US10068037B2/en
Publication of JP2016125934A publication Critical patent/JP2016125934A/ja
Application granted granted Critical
Publication of JP6458501B2 publication Critical patent/JP6458501B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)

Description

本発明は、シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置に関する。
数値計算を用いて流体である水や空気の流れを調べる流体解析や、圧縮された弾性体であるゴムの振る舞いを調べる弾性体解析などの、連続体の運動の解析において、粒子法と呼ばれる技術がある。連続体は、物質の巨視的な性質を保持した質点の集合体であり、質点は物質の存在する空間内に連続的に、密に分布している。粒子法は、連続体を粒子の分布で実現し、連続体の運動を粒子の運動として解析する解析手法である。粒子法としては、ＳＰＨ（Smoothed Particles Hydrodynamics）法やＭＰＳ（Moving Particles Semi-implicit）法が知られている。また、粒子法において粉流体の挙動を表現する手法としては、個別要素法と呼ばれる方法がある。
粒子分布を用いて解析対象を表現する場合、固定境界が設定される。固定境界を実現する手法には様々なものがある。
粒子法において固定境界を実現する標準的な方法としては、以下の３つの方法がある。第１の方法では、境界に沿って仮想的な流体粒子（以下、「境界粒子」と呼ぶ）を配置する。この境界粒子を固定境界に固着した流体とみなすことで、固定境界が実現される。第２の方法では、第１の方法と同様に配置された境界粒子と、解析対象である連続体を実現する粒子との間に、両粒子間の距離の関数で大きさが与えられ、両粒子の相対位置ベクトルに沿った方向の反発力が設定される。これにより、連続体を実現する粒子が固定境界の内側に運動の領域が限定される。第３の方法では、解析対象である連続体を実現する粒子に対して、予め設定された境界面を挟んで鏡映対称な位置に仮想的な連続体粒子を配置することで、固定境界が実現される。
固定境界を実現する他の方法としては、境界粒子を配置し、その境界粒子に対して一定以下の距離に近づいた連続体粒子に対して、境界面の法線に沿った向きの反発力を与えることで、固定境界を実現する方法がある。
また、固定境界をポリゴンの集合として表現する方法もある。この方法では、連続体粒子は、その重心座標が最も近いポリゴンとの間に、該連続体粒子と該ポリゴンを含む平面との距離に応じた反発力が与えられる。ただし、該連続体粒子と該ポリゴンを含む平面との距離が一定値以下になる場合にのみ、反発力が与えられる。
さらに、汎用性が高く、濡れ性モデルを統一的に取り扱うことができ、空間解像度が粗くても精度よく解析を行うことができる流体解析装置も考えられている。この流体解析装置は、流体を表す第１の粒子ごと、および流体に接する壁を表す第２の粒子ごとに、速度、位置、圧力を記録し、流体と壁との接触角θと、各粒子が他の粒子から作用力を受ける範囲に含まれる粒子を特定する。そして流体解析装置は、第１の粒子は第１の関数で、第２の粒子は第１の関数に（１＋ｃｏｓθ）／２を乗じて作用力を算出する。
さらに、次の処理をコンピュータに実行させる連続体運動解析プログラムが考えられている。連続体運動解析プログラムに基づいて、コンピュータは、連続体を粒子として表し、固定境界を任意の形状を有する微小な領域（微小面要素）の集合として表し、各連続体粒子の影響範囲に含まれる微小面要素から連続体粒子に対して微小面要素の法線方向に働く反発力を算出する。そしてコンピュータは、各反発力を足し合わせて連続体粒子に働く固定境界からの力を求める。
計算コストを抑えると共に固定壁の境界条件について任意の形状に対応可能であり、さらに計算精度の良い流動解析方法も考えられている。この流動解析方法では、壁境界の距離が所定以下と判断された第１粒子の各々に対して最も近い壁境界の曲率が算出される。そして、算出された曲率に応じた拡大率で、壁境界を挟んで対称に、壁を表す複数の第２粒子が、境界条件として配置される。
特開２００８−１１１６７５号公報 国際公開第２０１２／０２５９９５号 特開２０１０−２４３２９３号公報
G.R.Liu, M.B.Liu, "Smoothed Particle Hydrodynamics:A Meshfree Particle Method", World Scientific Pub Co Inc., 2003年10月, pp.138-141 M.G.Gesteira, B.D.Rogers, R.A.Dalrymplem, A.J.C.Crespo, M.Narayanaswamy, "User Guide for the SPHysics Code v1.4",2009年2月, pp.16-19 原田隆宏、越塚誠一、河口洋一郎、「流体と布のリアルタイム連成シミュレーション」、情報処理学会研究報告、グラフィクスとＣＡＤ研究会報告、社団法人情報処理学会、2007年11月12日、p.13-18
従来技術のうち、境界面を微小面要素の集合として表現する固定境界の表現方法は、他の表現方法に比べ、解析を行うために用いるメモリ容量および解析の計算時間を軽減することができ、適用場面に関する制約も少ないという利点がある。しかし、微小面要素から一定距離以内に近づいた連続体粒子に反発力を及ぼすため、一定距離よりも狭い隙間を有する領域内での連続体の運動解析を行うことができない。境界要素からの反発力が及ぶ距離は、例えば、数値計算上の空間解像度に合わせて平均粒子間隔程度とすることが多い。このとき、粒子間隔が１ｃｍ程度のときには、５ｍｍ程度の間隔で空いた隙間に粒子が入り込む過程を扱うことなどができない。
１つの側面では、本件は、狭い隙間を有する領域内での粒子法による解析を可能とすることを目的とする。
１つの案では、コンピュータに以下の処理を実行させるシミュレーションプログラムが提供される。コンピュータは、連続体を表す複数の粒子の位置と、第１の領域と第２の領域との境界を表す複数の要素の位置とに基づいて、複数の粒子それぞれが第２の領域内にあるか否かを判断する。次にコンピュータは、第２の領域内にある粒子に関して、複数の要素で表される境界と粒子との距離に基づいて、第１の領域側に向かう力を算出する。そしてコンピュータは、第２の領域内にある粒子に力を作用させて、複数の粒子の運動を解析する。
１態様によれば、狭い隙間を有する領域内での粒子法による解析が可能となる。
第１の実施の形態に係るシミュレーション装置の構成例を示す図である。 第２の本実施の形態に用いるコンピュータのハードウェアの一構成例を示す図である。 コンピュータのシミュレーション機能を示すブロック図である。 連続体粒子データの一例を示す図である。 円盤要素データの一例を示す図である。 交差円盤リストの一例を示す図である。 固定境界での反発力作用範囲の違いを示す図である。 凸に交差する円盤要素の例を示す図である。 凹に交差する円盤要素の例を示す図である。 シミュレーション手順の一例を示す第１のフローチャートである。 連続体粒子と円盤要素との位置関係の判断例を示す図である。 交差円盤リスト作成処理の手順の一例を示すフローチャートである。 ２つの円盤要素が共通部分を持つかどうかの判定例を示す図である。 シミュレーション手順の一例を示す第２のフローチャートである。 シミュレーション手順の一例を示す第３のフローチャートである。
以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔第１の実施の形態〕
第１の実施の形態は、流体や弾性体などの連続体を粒子法で解析するシミュレーション装置である。このシミュレーション装置は、連続体を表す粒子の運動解析時に、境界の裏側への粒子の侵入を許容し、境界の表側に近づいても反発させない。その代わり、シミュレーション装置は、境界の裏側に入った粒子に対して、表側に戻す反発力を加える。これにより、屈折する境界の狭い隙間にも粒子が入り込めるようになる。
図１は、第１の実施の形態に係るシミュレーション装置の構成例を示す図である。シミュレーション装置１０は、記憶部１１と演算部１２とを有する。
記憶部１１は、連続体を表す複数の粒子４〜７に関する粒子データ１１ａと、第１の領域１と第２の領域２との境界３を表す複数の境界要素３ａ〜３ｅに関する境界要素データ１１ｂとを記憶する。第１の領域１は、連続体で表される物質が存在する領域である。
第１の領域１と第２の領域２との境界３は、例えば、複数の境界要素３ａ〜３ｅを交差させながら組み合わせることで形成される。境界３は、要素間の交差部において屈折しており、狭い隙間が生じている場合もある。例えば図１の例では、境界要素３ｂと境界要素３ｄとの間が、狭い隙間となっている。
粒子データ１１ａには、例えば各粒子４〜７の位置を示す情報が含まれる。境界要素データ１１ｂには、各境界要素３ａ〜３ｅの位置を示す情報が含まれる。複数の境界要素３ａ〜３ｅは、例えば円盤形のような面状の要素である。複数の境界要素３ａ〜３ｅが円盤形の場合、境界要素データ１１ｂには、さらに各境界要素３ａ〜３ｅの大きさを示す情報（例えば半径または面積）や向きを示す情報（例えば法線ベクトル）が含まれる。
演算部１２は、複数の粒子４〜７の位置と、複数の境界要素３ａ〜３ｅの位置とに基づいて、複数の粒子４〜７それぞれが第２の領域２内にあるか否かを判断する。例えば複数の境界要素３ａ〜３ｅが円盤形であるとき、第１の領域１を向く面を表面とし、第２の領域２を向く面を裏面とする。このとき、演算部１２は、例えば、複数の粒子４〜７それぞれを注目粒子とする。そして演算部１２は、注目粒子の位置が裏面側となる境界要素が少なくとも１つ存在するかどうかを判断する。該当する境界要素が存在しない場合、演算部１２は、注目粒子が第１の領域１内にあると判断する。また該当する境界要素が存在する場合、演算部１２は、境界要素のうちの該注目粒子に最も近い最近接要素を特定する。そして演算部１２は、注目粒子の位置が、最近接要素に対して第１の領域１側に凸となるように交差する交差境界要素の表面側ではない場合、注目粒子が第２の領域２内にあると判断する。
第２の領域２内にある粒子４，７を検出すると、演算部１２は、それらの粒子４，７に関して、複数の境界要素３ａ〜３ｅで表される境界３と粒子４，７との距離に基づいて、第１の領域１側に向かう力４ａ，７ａを算出する。例えば演算部１２は、複数の境界要素３ａ〜３ｅが、第１の領域１を向く表面と第２の領域２を向く裏面とを有する円盤の場合、粒子４，７が裏面側にある境界要素のうちの、粒子に最も近い境界要素の表面側の法線方向の力を算出する。
また演算部１２は、粒子が、第１の領域１側に凹となるように交差する第１および第２の要素それぞれの裏面側にある場合、第１および第２の要素のうち、粒子との距離が長い方の要素における粒子との距離に基づいて、力を算出する。
力の算出後、演算部１２は、第２の領域２内にある粒子４，７に対して算出した力４ａ，７ａを作用させて、複数の粒子４〜７の運動を解析する。
このように、粒子が境界要素の裏側に回り込んだ時に反発力を与えることで、境界形状に細い隙間があっても、その隙間に容易に粒子が入り込むことができる。例えば粒子５は、いずれの要素の裏側でもないため、要素からの反発力を受けず、容易に隙間に入り込むことができる。その結果、流体や弾性体の粒子の流れ込みを扱う解析が可能となる。
また、複数の粒子４〜７は、第２の領域２側に入り込んだときに、例えば、境界３を構成する複数の境界要素３ａ〜３ｅのうち、粒子に対して裏側を見せているいずれかから反発する力を受けることとなる。例えば粒子４は、境界要素３ｄと境界要素３ｅとの裏面側にある。粒子４に最も近いのは境界要素３ｅであるため、粒子４は、境界要素３ｅから、例えば境界要素３ｅの法線方向の力４ａを受ける。このように、より近い境界要素３ｅから力４ａを受けることで、粒子４は、迅速に第１の領域１に戻すことができる。
なお、複数の境界要素３ａ〜３ｅを用いて境界３の形状を定義すると、要素の端が表現したい形状からはみ出ることがあり得る。第１の実施の形態では、要素同士の角度を用いて、反発力の与え方に修正を施すことで、要素の端において、不自然な反発力が発生することを防ぐことができる。例えば、粒子６の位置は、境界要素３ｄの裏側にあるが、その境界要素３ｄとの間で、第１の領域１側に凸に交差する境界要素３ｅの表側でもある。この場合、粒子６には反発力は働かない。すなわち、交差する要素の端部が飛び出していても、飛び出した部分の影響で粒子６に不自然な力が働くことが抑止される。なお粒子７の位置は、境界要素３ｃの裏側にあり、境界要素３ｃと交差する境界要素３ｂの表側にあるが、境界要素３ｃと境界要素３ｄとは、第１の領域１に対して凹に交差している。そのため、粒子７は、例えば境界要素３ｃの法線方向の力７ａを境界要素３ｃから受けることとなる。
しかも、要素間の交差部分において要素の端が飛び出ることが許容されることにより、複数の境界要素３ａ〜３ｅを用いて境界３の形状を表現する際に、要素間の接続情報を取り扱わずにすむ。その結果、形状を生成するときに煩雑な処理が減らされ、表面を隙間なく覆う境界３の作成が容易となる。
図１に示す第１の実施の形態に示すシミュレーション装置１０による連続体の運動解析における境界面の設定手法は、他の手法に比べ、さまざまな点で優れている。
例えば、境界粒子を配置し、その境界粒子に対して一定以下の距離に近づいた連続体粒子に対して、境界面の法線に沿った向きの反発力を与える手法では、滑らかな形状の境界を表現するためには多数の粒子を用いることとなる。粒子の量が多ければ、メモリ容量、計算時間のいずれも多く消費される。それに対して第１の実施の形態では、粒子で境界を構成する場合に使用する粒子よりも少ない量の要素で境界を設定でき、消費するメモリ量の削減、および計算時間の削減が可能となる。
しかも多数の粒子で境界を設定する手法の中には、境界により粒子に与える反発力は、例えば、連続体を構成する粒子の位置からｘ軸，ｙ軸それぞれの正・負の方向に１つずれた位置にある４つの境界粒子を用いて得られる法線ベクトル方向に設定されるものがある。この場合、境界面上の粒子分布は正則（隣接する粒子数が一定となること）な格子点、もしくは、それを連続的に変形させた位置に配置しなくてはならないという制約がある。その結果、例えば、テトラメッシュの重心位置に境界粒子を配置するなどした場合には応用することができないという問題がある。それに対して、第１の実施の形態では、境界の設定に関し、そのような制約はなく、柔軟性が高い。
また、固定境界をポリゴンの集合として表現する手法もあるが、この手法ではポリゴンのサイズが一様であると仮定している。そのため、ポリゴンのサイズが異なる状況では、ポリゴンの境界付近において不合理な反発力を生じうるという問題がある。それに対して、第１の実施の形態のシミュレーション装置１０では、境界を構成する要素のサイズが一定でなくても適切な反発力を生じさせることができる。
さらに、境界面に対して鏡映対称な位置に仮想的な連続体粒子を配置する手法がある。この手法では、鏡映対称な位置を計算するためには境界面を解析的に表現しなくてはならない。そのため、例えば、人手によって作成された、境界面の形状、および、その法線方向を表現する数式が特に与えられていないような境界形状に対しては適用できないという問題がある。それに対して、第１の実施の形態では、境界面の法線方向は、円盤形の要素データから一意に決定するため、数式によって与えずにすむ。
粒子に最も近い壁境界の曲率に応じた拡大率で、壁境界を挟んで対称に、壁を表す複数の第２粒子を配置する手法では、境界面が動くとき、距離関数の算出を毎ステップ行うこととなり、計算量が過大となる。そのため、数値計算に要する時間が増大するという問題がある。それに対して、第１の実施の形態では、境界面が動く場合、境界データを更新すればよく、数値計算に要する時間の増大は微小である。
固定境界を任意の形状を有する微小な領域（微小面要素）の集合として表す従来の手法では、境界要素から一定距離以内に近づいた粒子に反発力を及ぼす。そのため、その距離よりも狭い隙間を扱うことができない。境界要素からの反発力が及ぶ距離は数値計算上の空間解像度に合わせて平均粒子間隔程度とすることが多いが、例えば、粒子間隔が１ｃｍ程度のときには、５ｍｍ程度の間隔で空いた隙間に粒子が入り込む過程を扱うことなどができない。それに対して、第１の実施の形態では、狭い隙間にも粒子を入り込ませることができ、狭い隙間を有する領域における連続体の運動を、適切に解析することができる。
なお、演算部１２は、例えばシミュレーション装置１０が有するプロセッサにより実現することができる。また、記憶部１１は、例えばシミュレーション装置１０が有するメモリにより実現することができる。
〔第２の実施の形態〕
次に第２の実施の形態について説明する。第２の実施の形態は、粒子法を用いた数値計算により、流体または弾性体の運動解析シミュレーションを行うものである。なお第２の実施の形態では、円形の円盤要素で境界面が設定される。
図２は、第２の本実施の形態に用いるコンピュータのハードウェアの一構成例を示す図である。コンピュータ１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。
メモリ１０２は、コンピュータ１００の主記憶装置として使用される。メモリ１０２には、プロセッサ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、メモリ１０２には、プロセッサ１０１による処理に必要な各種データが格納される。メモリ１０２としては、例えばＲＡＭ（Random Access Memory）などの揮発性の半導体記憶装置が使用される。
バス１０９に接続されている周辺機器としては、ＨＤＤ（Hard Disk Drive）１０３、グラフィック処理装置１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７およびネットワークインタフェース１０８がある。
ＨＤＤ１０３は、内蔵したディスクに対して、磁気的にデータの書き込みおよび読み出しを行う。ＨＤＤ１０３は、コンピュータ１００の補助記憶装置として使用される。ＨＤＤ１０３には、ＯＳのプログラム、アプリケーションプログラム、および各種データが格納される。なお、補助記憶装置としては、フラッシュメモリなどの不揮発性の半導体記憶装置を使用することもできる。
グラフィック処理装置１０４には、モニタ２１が接続されている。グラフィック処理装置１０４は、プロセッサ１０１からの命令に従って、画像をモニタ２１の画面に表示させる。モニタ２１としては、ＣＲＴ（Cathode Ray Tube）を用いた表示装置や液晶表示装置などがある。
入力インタフェース１０５には、キーボード２２とマウス２３とが接続されている。入力インタフェース１０５は、キーボード２２やマウス２３から送られてくる信号をプロセッサ１０１に送信する。なお、マウス２３は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。
光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク２４に記録されたデータの読み取りを行う。光ディスク２４は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク２４には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）などがある。
機器接続インタフェース１０７は、コンピュータ１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置２５やメモリリーダライタ２６を接続することができる。メモリ装置２５は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ２６は、メモリカード２７へのデータの書き込み、またはメモリカード２７からのデータの読み出しを行う装置である。メモリカード２７は、カード型の記録媒体である。
ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。
以上のようなハードウェア構成によって、第２の実施の形態の処理機能を実現することができる。なお、第１の実施の形態に示したシミュレーション装置１０も、図２に示したコンピュータ１００と同様のハードウェアにより実現することができる。シミュレーション装置１０内の要素のうち、演算部１２は、コンピュータ１００のプロセッサ１０１によって実現できる。またシミュレーション装置１０内の要素のうち、記憶部１１は、コンピュータ１００のメモリ１０２またはＨＤＤ１０３によって実現できる。
コンピュータ１００は、例えばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態の処理機能を実現する。コンピュータ１００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。例えば、コンピュータ１００に実行させるプログラムをＨＤＤ１０３に格納しておくことができる。プロセッサ１０１は、ＨＤＤ１０３内のプログラムの少なくとも一部をメモリ１０２にロードし、プログラムを実行する。またコンピュータ１００に実行させるプログラムを、光ディスク２４、メモリ装置２５、メモリカード２７などの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、例えばプロセッサ１０１からの制御により、ＨＤＤ１０３にインストールされた後、実行可能となる。またプロセッサ１０１が、可搬型記録媒体から直接プログラムを読み出して実行することもできる。
図３は、コンピュータのシミュレーション機能を示すブロック図である。コンピュータ１００は、記憶部１１０、シミュレーション制御部１２０、運動計算部１３０、および反発力計算部１４０を有する。
記憶部１１０は、シミュレーションに用いるデータを記憶する。例えば記憶部１１０は、連続体粒子データ群１１１、円盤要素データ群１１２、および交差円盤リスト１１３を記憶する。記憶部１１０としては、例えばメモリ１０２またはＨＤＤ１０３の記憶領域の一部が使用される。
連続体粒子データ群１１１は、連続体粒子ごとに用意された単位データ（連続体粒子データ）の集合である。また、連続体粒子データは流体や弾性体をモデル化したものであり、流体・弾性体などの連続体を表す各粒子の状態を示すデータである。連続体粒子データ群１１１は、シミュレーション制御部１２０によりシミュレーション開始時に初期値が設定され、シミュレーションの進行に伴い更新される。
円盤要素データ群１１２は、円盤要素ごとに用意された単位データ（円盤要素データ）の集合であり、流体を格納した容器の壁などを微小円盤の集合として境界の壁全体をモデル化したものである。円盤要素データは、境界面を形成する各円盤要素の形状を示すデータである。円盤要素データ群１１２は、シミュレーション制御部１２０によりシミュレーション開始時に設定される。
交差円盤リスト１１３は、ある円盤要素と交差する他の円盤要素のリストである。交差円盤リスト１１３は、反発力計算部１４０による反発力計算時に作成される中間データである。
シミュレーション制御部１２０は、運動計算部１３０および反発力計算部１４０と連携して、連続体の運動解析のシミュレーションを行う。例えばシミュレーション制御部１２０は、シミュレーションの時間ステップごとに、その時間ステップにおける連続体粒子の状態を計算する。そしてシミュレーション制御部１２０は、連続体粒子の運動（速度や加速度）については、運動計算部１３０に計算させる。またシミュレーション制御部１２０は、連続体粒子に対して円盤要素から働く反発力については、反発力計算部１４０に計算させる。
運動計算部１３０は、流体または弾性体に働く物理法則に従って、シミュレーション上のタイムステップごとの連続体粒子の位置や物理量を計算する。運動計算部１３０は、計算結果に基づいて、連続体粒子データ群１１１内の各連続体粒子データを更新する。
反発力計算部１４０は、連続体粒子の運動に対して固定境界条件を適用するため、円盤要素に一定以上近づいた連続体粒子に対して、円盤要素の法線方向に加える反発力を計算する。例えば反発力計算部１４０は、連続体粒子が円盤要素の裏側にあるときに、その円盤要素の表側に向かう反発力を、その連続体粒子に働かせる。
なお、図３に示した各要素の機能は、例えば、その要素に対応するプログラムモジュールをコンピュータ１００に実行させることで実現することができる。
次に、記憶部１１０内のデータについて詳細に説明する。
図４は、連続体粒子データの一例を示す図である。連続体粒子データ群１１１に含まれる連続体粒子データ１１１ａ，１１１ｂ，１１１ｃ，・・・には、例えば連続体粒子番号、中心座標データ、速度データ、加速度データ、およびその他のデータが含まれる。連続体粒子番号は、連続体粒子の識別番号である。中心座標データは、連続体粒子の中心（例えば重心）位置の座標である。速度データは、連続体粒子の移動速度である。加速度データは、連続体粒子の加速度である。その他のデータは、例えば連続体粒子で表現されている流体の物性値などである。例えば、影響半径、密度、質量、変形勾配テンソル、材料物性、温度などのデータが、その他のデータに含まれる。
図５は、円盤要素データの一例を示す図である。円盤要素データ群１１２に含まれる円盤要素データ１１２ａ，１１２ｂ，１１２ｃ，・・・には、例えば円盤要素番号、中心座標データ、法線ベクトルデータ、面積データなどが含まれる。円盤要素番号は、円盤要素の識別番号である。中心座標データは、円盤要素の中心（例えば円盤要素の重心）位置の座標である。法線ベクトルデータは、円盤要素の表側に向かう大きさ１の法線ベクトルを示すベクトルデータである。面積データは、円盤要素の面積を示すデータである。面積データに基づいて、円盤要素の半径を算出できる。
図６は、交差円盤リストの一例を示す図である。交差円盤リスト１１３には、判定対象となっている円盤要素に交差する他の円盤要素それぞれに関するエントリ１１３ａ，１１３ｂ，・・・が登録されている。各エントリには、円盤要素番号と角度フラグデータとが含まれる。円盤要素番号は、判定対象の円盤要素に交差する他の円盤要素の識別番号である。角度フラグデータは、対応する円盤要素番号で示される円盤要素と、判定対象の円盤要素（最近接要素）との交差する位置における表側の角度が、凸なのか凹なのかを示すデータである。交差する位置の表側の角度が１８０度以上であれば、凸に交わっており、交差する位置の表側の角度が１８０度未満であれば、凹に交わっている。
このような機能を有するコンピュータ１００により、連続体の運動解析シミュレーションが行われる。なお、第２の実施の形態におけるシミュレーションでは、連続体粒子が固定境界の裏側に入り込むことを許容している。そして、固定境界の境界面が円盤要素で表され、円盤要素の裏側に入り込んだ連続体粒子に対して、反発力が加えられる。これにより、固定境界の表側で反発力を生じさせた場合と異なり、極めて狭い空間にも連続体粒子が侵入できるようになる。
図７は、固定境界での反発力作用範囲の違いを示す図である。図７では、上段に、固定境界の表側の連続体粒子に反発力を加えた例を示し、下段に、固定境界の裏側の連続体粒子に反発力を加えた例を示している。固定境界は、複数の円盤要素３１〜３５によって境界面が表現されている。図７の例では、円盤要素３１〜３５で表された境界面の上が表側、下が裏側である。
固定境界の表側の連続体粒子に反発力を加える場合、円盤要素３１〜３５の表側の所定の距離ｄ以下の範囲が、反発力作用範囲４１となる。反発力作用範囲４１外にある連続体粒子５１には反発力は加えられないが、反発力作用範囲４１内にある連続体粒子５２には、円盤要素３５から遠ざかる方向への反発力４２が加えられる。
ここで、反発力作用範囲４１を定める距離ｄに比べて狭い幅（例えば距離ｄの２倍以下）の領域４３があると、その領域４３内に連続体粒子５１，５２が入り込むことは困難となる。図７の例では、円盤要素３２〜３４に囲まれた領域が、狭い領域４３である。また、連続体粒子５１，５２が狭い領域４３に入り込んだとしても、常にいずれかの円盤要素から反発力を受けることとなり、運動を正しく解析することができない。
他方、固定境界の裏側の連続体粒子に反発力を加える場合、円盤要素３１〜３５の裏側の所定の距離ｄ以下の範囲が、反発力作用範囲４４となる。反発力作用範囲４４外にある連続体粒子５３には反発力は加えられないが、反発力作用範囲４４内にある連続体粒子５４には、円盤要素３５の表側に向かう反発力４５が加えられる。
この場合、幅の狭い領域４３は、反発力作用範囲４４外である。そのため、連続体粒子５３，５４は、領域４３内に容易に入り込むことができる。しかも、狭い領域４３に入り込んだ連続体粒子５３に対して、反発力が作用しないため、運動を正しく解析することができる。
なお、固定境界を形成する複数の円盤要素は、隣接する円盤要素同士が交差する。そのため、各円盤要素には、隣接する円盤要素と交差する部分からはみ出す部分がある。三次元空間内で円盤要素を組み合わせた際に生じる円盤要素の端のはみ出しは、凸に交差する場合と凹に交差する場合とに分けて、以下の方針で処理される。
図８は、凸に交差する円盤要素の例を示す図である。図８の例では、円盤要素３６と円盤要素３７とが凸に交差している。円盤要素３６の中心位置は、位置ベクトルｒj1で表される。円盤要素３６の法線ベクトルはｎj1である。同様に円盤要素３７の中心位置は、位置ベクトルｒj2で表される。円盤要素３７の法線ベクトルはｎj2である。円盤要素３７の位置ベクトルｒj2から、円盤要素３６の位置ベクトルｒj1を減算することで、円盤要素３６から円盤要素３７への相対位置ベクトルｒj1j2が得られる。
図８に示したように円盤要素３６，３７が凸に交差する状況では、各円盤要素３６，３７の裏側に入り込んだ連続体粒子ｒi1は、最も鉛直距離の近い円盤要素３６から反発力を受ける。ただし、最も鉛直距離の近い円盤要素３６の裏側ではあるが、円盤要素３６に凸に交わる円盤要素３７の表側にある連続体粒子ｒi2は、いずれの円盤要素３６，３７からも反発力を受けないこととする。
このように、連続体粒子が円盤要素の裏側にある場合でも、隣接する他の円盤要素の表側となる場合には、その連続体粒子には反発力を加えないことで、本来の固定境界の内側（円盤要素の表面側）にある連続体粒子には、円盤要素からの反発力が働かなくなる。その結果、交差する円盤要素における交差位置よりも先の部分においても、流体または弾性体の運動を正しく解析できる。
２つの円盤要素３６，３７の交わる角度の凹凸の判定は、２つの円盤要素３６，３７間の相対位置ベクトルｒj1j2と、法線ベクトルｎj1，ｎj2それぞれとの内積の符号を用いて判断できる。
例えば反発力計算部１４０は、円盤要素３６，３７それぞれの円盤要素データから、法線ベクトルｎj1，ｎj2を示す法線ベクトルデータを取得する。また反発力計算部１４０は、円盤要素３６，３７それぞれの円盤要素データから面積データを取得し、その面積データから円盤要素３６，３７の半径Ｒj1，Ｒj2を求める。さらに反発力計算部１４０は、円盤要素３６から円盤要素３７それぞれの中心座標データの各座標成分の値の差を取って、相対位置ベクトルｒj1j2とする。
ここで、次の条件の両方を満たすとき、円盤要素３６と円盤要素３７とは、凸なる角度を持っていると判断できる。
条件１：ｒj1j2・ｎj1＜０
条件２：−ｒj1j2・ｎj2＜０
また、いずれの条件も満たさないとき、円盤要素３６と円盤要素３７とは凹なる角度を持っていると判断できる。交差する２円盤要素が条件の片方のみを満たし、もう片方を満たさないとき、これらの２要素によって境界面の表裏を整合的に表現することはできないので、そのような場合は、処理を止めるか、凹凸いずれかの場合に準じた処理を行う。
図９は、凹に交差する円盤要素の例を示す図である。図９の例では、円盤要素３６と円盤要素３７とが凹に交差している。この状況において、各円盤要素３６，３７それぞれ裏側に入り込んだ連続体粒子ｒi3は、鉛直距離の遠い方から反発力を受けることとする。例えば連続体粒子ｒi3から円盤要素３６への鉛直距離がｒn1であり、連続体粒子ｒi3から円盤要素３７への鉛直距離がｒn2であるとする。図９の例では、ｒn1＜ｒn2である。そのため、連続体粒子ｒi3は、円盤要素３７から、円盤要素３７の法線ベクトルｎj2の正の方向への反発力を受ける。
このように、連続体粒子が複数の円盤要素の裏側にある場合、鉛直距離が遠い方の円盤要素から反発力を働かせることで、連続体粒子が固定境界を超えて入り込む深さを、迅速に減少させることができる。
次に、シミュレーションの手順について詳細に説明する。
図１０は、シミュレーション手順の一例を示す第１のフローチャートである。
［ステップＳ１０１］シミュレーション制御部１２０は、入力データを取得する。入力データには、数値計算を行う直接の対象である連続体粒子データや、連続体粒子の運動の境界条件を設定する円盤要素データが含まれる。入力データは、例えば、ユーザ操作によってコンピュータ１００に入力される。シミュレーション制御部１２０は、入力データを、記憶部１１０に格納する（図４、図５参照）。
シミュレーション制御部１２０は、入力データの取得後、シミュレーションの時間ステップごとに、以下の処理を実行する。
［ステップＳ１０２］運動計算部１３０は、連続体粒子に設定された物理モデルに従って、連続体粒子の運動を解く。例えば、運動計算部１３０は、流体をモデル化した連続体粒子に対しては、ナヴィエ・ストークス方程式や連続方程式を解いて、加速度や密度の時間微分を求める。
［ステップＳ１０３］反発力計算部１４０は、反発力を計算していない連続体粒子の中から、反発力の計算対象とする連続体粒子を１つ選択する。以下、選択した連続体粒子を、注目連続体粒子と呼ぶ。
［ステップＳ１０４］反発力計算部１４０は、注目連続体粒子に対して、以下の条件を満たした円盤要素を同定し、最近接円盤とする。
条件１：円盤要素から見て、注目連続体粒子の位置が法線の方向（表側）の逆側（裏側）に存在する。
条件２：円盤要素と法線方向の作る円筒の内部に存在する。
条件３：上記２つの条件を共に満たす円盤要素のなかで、注目連続体粒子との鉛直距離が最も小さい。
注目連続体粒子と円盤要素が上記条件を満たしているかどうかの判断は、次の通り行う。
図１１は、連続体粒子と円盤要素との位置関係の判断例を示す図である。ここで、注目連続体粒子６１の番号をｉ（ｉは１以上の整数），円盤要素６２の番号をｊ（ｊは１以上の整数）とする。反発力計算部１４０は、注目連続体粒子６１の位置から、円盤要素６２の中心への相対位置ベクトルを求める。相対位置ベクトルは、注目連続体粒子６１の中心座標データの各座標成分と円盤要素６２の座標データの各座標成分との差を取ることで求めることができる。反発力計算部１４０は、求められた相対位置ベクトルを、ｒijとする。
次にｊ番の円盤要素６２の法線ベクトルデータの値をｎj, ｊ番の円盤要素６２の面積データの値をＡjとする。反発力計算部１４０は、ｒijのｎj方向成分ｒnijとｒijのｎjに垂直な成分ｒtijとを次の式から求める。
ｒnij＝（ｒij・ｎj）ｎj
ｒtij＝ｒij−ｒnij
反発力計算部１４０は、円盤要素６２の半径ＲjをＲj＝（Ａj／π）1/2として求める（πは円周率）。ここで、「条件１」が満たされているかどうかを、ｒij・ｎj＜０であるかどうかによって判断することができる。すなわち、ｒij・ｎj＜０であれば、注目連続体粒子６１は円盤要素６２の表側にあり、「条件１」は満たされない。ｒij・ｎj＜０でなければ、注目連続体粒子６１は円盤要素６２の裏側にあり、「条件１」は満たされる。なお、この例では、ｒij・ｎj＝０の場合も、注目連続体粒子６１が円盤要素６２の裏側にあるとみなしている。
「条件２」が満たされているかどうかは、Ｒj＞｜ｒtij｜という条件式が満たされているかどうかで判断できる。すなわち、Ｒj＞｜ｒtij｜が満たされれば、「条件２」が満たされる。この例では、Ｒj＝｜ｒtij｜の場合には、「条件２」が満たされないこととなるが、Ｒj＝｜ｒtij｜の場合、「条件２」が満たされるものとしてもよい。
ｉ番目の注目連続体粒子６１に対して、「条件１」と「条件２」との２つの条件を共に満たす円盤要素のなかで、｜ｒnij｜が最も小さい円盤要素が、「条件３」を満たす最近接円盤である。
このような最近接円盤の判断が、すべての連続体粒子に対して行われる。なお、すべての円盤要素の表側にある連続体粒子に対しては、最近接円盤は存在しない。少なくとも１つの円盤要素の裏側に存在する連続体粒子に対しては、最近接円盤が求められる。
以下、図１０の説明に戻る。
［ステップＳ１０５］反発力計算部１４０は、注目連続体粒子に対する最近接円盤が存在するか否かを判断する。最近接円盤が存在しない場合、反発力計算部１４０は、注目連続体粒子がいずれの円盤要素からも反発力を受けないと判断し、処理をステップＳ１２１（図１５参照）に進める。最近接円盤が存在すれば、処理がステップＳ１０６に進められる。
［ステップＳ１０６］反発力計算部１４０は、最近接円盤に交差する円盤要素のリストである交差円盤リスト１１３（図６参照）を作成する。交差円盤リスト作成処理の詳細は後述する（図１２参照）。その後、処理がステップＳ１１１（図１４参照）に進められる。
次に、交差円盤リスト１１３作成処理について詳細に説明する。
図１２は、交差円盤リスト作成処理の手順の一例を示すフローチャートである。
［ステップＳ２０１］反発力計算部１４０は、最近接円盤を円盤要素ｊ1とする。
［ステップＳ２０２］反発力計算部１４０は、円盤要素データ群１１２に含まれる円盤要素データに示される円盤要素を１つ選択する。
［ステップＳ２０３］反発力計算部１４０は、選択した円盤要素を、円盤要素ｊ2とする。
［ステップＳ２０４］反発力計算部１４０は、円盤要素ｊ1を含む平面Ｐ1と円盤要素ｊ2とが共通部分を持つかどうかを判断する。
図１３は、２つの円盤要素が共通部分を持つかどうかの判定例を示す図である。これは｜ｒj1j2・ｎj1｜＜Ｒ12ｓｉｎθ12であれば、円盤要素ｊ1を含む平面Ｐ1と円盤要素ｊ2とが共通部分を持つと判断できる。ここでθ12は２つの円盤要素それぞれの法線ベクトルｎj1，ｎj2がなす角であり、θ12＝ａｒｃｃｏｓｎj1・ｎj2として計算できる。｜ｒj1j2・ｎj1｜は、円盤要素ｊ2の中心から平面Ｐ1までの鉛直距離を示している。Ｒ12ｓｉｎθ12は、円盤要素ｊ2の外縁が、円盤要素ｊ1の法線方向に、円盤要素ｊ2の中心からどの程度の距離まで届くかを示している。
以下、図１２の説明に戻る。
円盤要素ｊ1を含む平面Ｐ1と円盤要素ｊ2とが共通部分を持つ場合、処理がステップＳ２０５に進められる。また、共通部分を持たない場合、円盤要素は互いに交わらないと判断され、選択した円盤要素を交差円盤リスト１１３に加えずに、処理がステップＳ２１１に進められる。
［ステップＳ２０５］反発力計算部１４０は、円盤要素ｊ2を含む平面Ｐ2と円盤要素ｊ1とが共通部分を持つかどうかを判断する。共通部分を持つ場合、処理がステップＳ２０６に進められる。また、共通部分を持たない場合、円盤要素は互いに交わらないと判断され、選択した円盤要素を交差円盤リスト１１３に加えずに、処理がステップＳ２１１に進められる。
［ステップＳ２０６］反発力計算部１４０は、ステップＳ２０４，Ｓ２０５のいずれの場合においても共通部分を持つなら、円盤要素ｊ2の外縁と平面Ｐ1との交点（図１３参照）を求める。
これは未知ベクトルｒに関する連立方程式を解くことで求めることができる。ここで、円盤要素の中心座標データからｒj1，ｒj2、法線ベクトルデータからｎj1，ｎj2、面積データからＲj2が求まり、未知数としてｒを含む、以下の方程式をたてることができる。
（ｒ−ｒj1）・ｎj1＝０
（ｒ−ｒj2）・ｎj2＝０
（ｒ−ｒj2）2＝Ｒj2 2
これを解くと多くとも２つの交点座標が求まる。
［ステップＳ２０７］反発力計算部１４０は、求められた交点のうちの少なくとも１つが円盤要素ｊ1に含まれるかどうかを判断する。例えば、円盤要素ｊ1の中心から交点までの距離が、円盤要素ｊ1の半径以下であれば、その交点は円盤要素ｊ1に含まれる。少なくとも１つの交点が円盤要素ｊ1に含まれる場合、円盤要素ｊ1と円盤要素ｊ2とは交わるものと判断され、処理がステップＳ２１０に進められる。いずれの交点も円盤要素ｊ1に含まれない場合、処理がステップＳ２０８に進められる。
［ステップＳ２０８］反発力計算部１４０は、円盤要素ｊ１ の外縁と平面Ｐ２との交点を求める。
［ステップＳ２０９］反発力計算部１４０は、求められた交点のうちの少なくとも１つが円盤要素ｊ2に含まれるかどうかを判断する。少なくとも１つの交点が円盤要素ｊ2に含まれる場合、円盤要素ｊ1と円盤要素ｊ2とは交わるものと判断され、処理がステップＳ２１０に進められる。いずれの交点も円盤要素ｊ2に含まれない場合、円盤要素ｊ1と円盤要素ｊ2とは交わらないと判断され、処理がステップＳ２１１に進められる。
［ステップＳ２１０］反発力計算部１４０は、交差すると判定されたときの円盤要素ｊ2を、交差円盤リスト１１３に登録する。例えば反発力計算部１４０は、円盤要素ｊ2の円盤要素番号と、円盤要素ｊ1と円盤要素ｊ2とが凸に交わるのか凹に交わるのかを示す角度フラグデータとの組を、交差円盤リスト１１３に登録する。なお、２つの円盤要素ｊ1，ｊ2の交わる角度の凹凸の判定は、２つの円盤要素ｊ1，ｊ2間の相対位置ベクトルｒj1j2と、法線ベクトルｎj1，ｎj2それぞれとの内積の符号を用いて判断できる。
［ステップＳ２１１］反発力計算部１４０は、円盤要素データ群１１２内に、未選択の円盤要素があるか否かを判断する。未選択の円盤要素があれば、処理がステップＳ２０２に進められる。未選択の円盤要素がなければ、交差円盤リスト作成処理が終了する。
このようにして、最近接円盤と交差している円盤要素の番号と、その円盤要素と最近接円盤との交わる角を示すフラグ変数が登録された交差円盤リスト１１３が作成される。作成された交差円盤リスト１１３を用いて、注目連続体要素に対する反発力が計算される。
図１４は、シミュレーション手順の一例を示す第２のフローチャートである。図１４に示す処理では、注目連続体粒子に対して反発力を与える円盤要素（反発力対象要素）が決定される。
［ステップＳ１１１］反発力計算部１４０は、注目連続体粒子に対する最近接円盤を、反発力仮対象要素とする。以下の説明では、ｉ番目の連続体粒子が注目連続体粒子であるものとし、注目連続体粒子ｉと表す。
例えば、ｉ番目の連続体粒子の最近接円盤を、番号ｊ0で表し、この最近接円盤ｊ0に対して作成した交差円盤リストに属する円盤要素をｊ1〜ｊnとする。反発力計算部１４０は、まず、ｉ番目の連続体粒子の最近接円盤ｊ0を反発力仮対象要素とし、その番号ｊ0と連続体粒子との鉛直距離を、メモリに記憶する。
［ステップＳ１１２］反発力計算部１４０は、最近接円盤ｊ0に交差する未選択の円盤要素ｊ1〜ｊnを、交差円盤リスト１１３から１つ選択する。以下、選択した円盤要素を、注目円盤要素ｊkとする。
［ステップＳ１１３］反発力計算部１４０は、注目円盤要素ｊkが、最近接円盤ｊ0と凸に交わり、かつ、注目連続体粒子ｉが注目円盤要素ｊkの表側に存在するという条件が満たされるかどうかを調べる。注目円盤要素ｊkが最近接円盤ｊ0と凸に交わるかどうかは、交差円盤リスト１１３の角度フラグデータを用いて判断できる。また注目連続体粒子ｉが注目円盤要素ｊkの表側に存在するかどうかは、注目連続体粒子ｉから注目円盤要素ｊkへの相対位置ベクトルｒijkと注目円盤要素ｊkの法線ベクトルｎjkの内積の符号に基づいて判断できる。内積の符号が負であれば、注目連続体粒子ｉが注目円盤要素ｊkの表側にある。条件が満たされた場合、注目連続体粒子ｉは注目円盤要素ｊ k から反発力を受けないと判断され、処理がステップＳ１２１（図１５参照）に進められる。条件が満たされない場合、処理がステップＳ１１４に進められる。
［ステップＳ１１４］反発力計算部１４０は、注目円盤要素ｊkが最近接円盤ｊ0と凹に交わり、かつ注目連続体粒子ｉが注目円盤要素ｊkの法線の裏側に存在するという条件が満たされるかどうかを調べる。注目円盤要素ｊkが最近接円盤ｊ0と凹に交わるかどうかは、交差円盤リスト１１３の角度フラグデータを用いて判断できる。注目連続体粒子ｉが注目円盤要素ｊkの法線の裏側に存在するかどうかは、注目連続体粒子ｉから注目円盤要素ｊkへの相対位置ベクトルｒijkと円盤要素jkの法線ベクトルｎjkの内積の符号に基づいて判断できる。内積の符号が正であれば、注目連続体粒子ｉが注目円盤要素ｊkの裏側にある。条件が満たされた場合、処理がステップＳ１１５に進められる。条件が満たされない場合、処理がステップＳ１１７に進められる。
［ステップＳ１１５］反発力計算部１４０は、注目連続体粒子ｉと注目円盤要素ｊkとの鉛直距離が、反発力仮対象要素と注目連続体粒子ｉとの鉛直距離よりも大きいかどうかを調べる。注目連続体粒子ｉと注目円盤要素ｊkと鉛直距離は、注目連続体粒子ｉから注目円盤要素ｊkへの相対位置ベクトルｒijkと、注目円盤要素ｊkの法線ベクトルｎjkとの内積の絶対値｜ｒijk・ｎjk｜である。内積の絶対値と反発力仮対象要素との鉛直距離を比較することで、どちらの鉛直距離が大きいのかを判断できる。注目連続体粒子ｉとの鉛直距離が、注目円盤要素ｊkの方が大きければ、処理がステップＳ１１６に進められる。そうでなければ、処理がステップＳ１１７に進められる。
［ステップＳ１１６］反発力計算部１４０は、現在の注目円盤要素ｊkを、反発力仮対象要素とする。そして反発力計算部１４０は、注目円盤要素ｊkと注目連続体粒子ｉとの鉛直距離を、メモリに記憶しておく。
［ステップＳ１１７］反発力計算部１４０は、交差円盤リスト１１３内に、未選択の円盤要素があるか否かを判断する。未選択の円盤要素があれば、処理がステップＳ１１２に進められる。すべての円盤要素が選択済みであれば、処理がステップＳ１１８に進められる。
［ステップＳ１１８］反発力計算部１４０は、交差円盤リスト１１３内の全円盤要素に対してＳ１１３〜Ｓ１１６の処理を適用し終えた時点で反発力仮対象要素となっている円盤要素を、反発力対象要素とする。そして反発力計算部１４０は、反発力対象要素である円盤要素からの反発力を算出する。
ｉ番目の連続体粒子とｊ番目の円盤要素の間に働く力は、円盤要素の法線に沿って測った鉛直距離｜ｒnij｜に基づいて算出される。反発力の働く方向は、円盤要素の法線方向である。例えば、反発力ｆijは、以下の式で表される。
ｆij＝ｎjｆ0,i｛ｑ／（１−ｑ）｝
ここで、ｑ＝｜ｒnij｜／rcutoffである。rcutoffは円盤表面からその距離以上は裏側に入り込んでいかないという距離スケールを表す長さであり、予め設定される定数である。rcutoffとしては、例えば連続体粒子の持つ空間解像度程度の値が用いられる。空間解像度は、流体または弾性体の運動を、空間的にどの程度細かく表せるかを示す値である。空間解像度は、例えば流体または弾性体を構成する連続体粒子間の平均距離である。
ｆ0,iは、基準となる反発力の大きさであり、予め設定される定数である。ｆ0,iとしては、例えばｍiｃs,i 2／rcutoffとすればよい。ただし、ｍiは連続体粒子ｉの質量であり、ｃs,iは連続体粒子ｉの音速である。こうすることで、極めて大きな運動エネルギーを持った連続体粒子であっても、固定境界の裏側の反発力作用範囲をつきぬけるという事態が起こらないようにできる。
［ステップＳ１１９］反発力計算部１４０は、求めた反発力ｆijから定まる加速度を、ステップＳ１０２において既に求めてあった注目連続体粒子の加速度に足し合わせる。その後、処理がステップＳ１２１（図１５参照）に進められる。
図１５は、シミュレーション手順の一例を示す第３のフローチャートである。
［ステップＳ１２１］反発力計算部１４０は、未選択の連続体粒子があるか否かを判断する。未選択の連続体粒子があれば、処理がステップＳ１０３（図１０参照）に進められる。すべての連続体粒子が選択され、反発力の計算が完了していれば、処理がステップＳ１２２に進められる。
［ステップＳ１２２］運動計算部１３０は、反発力に基づいて求められた加速度を含む時間微分項を用いて、すべての連続体粒子の物理量（速度、加速度など）を時間積分し、１タイムステップ分の時間進行後の連続体粒子の位置や物理量を算出する。その後、シミュレーション制御部１２０が、シミュレーション上の時刻を時間刻みの分だけ進める。これにより、１タイムステップ分の計算が終了する。
［ステップＳ１２３］シミュレーション制御部１２０は、１タイムステップ分の計算結果データを外部ファイルとして出力する。出力されたファイルは、例えばＨＤＤ１０３に格納される。該出力ファイルは、別に用意されたプログラムで必要な数値や図を画面に表示したり、紙に印刷したりするときに利用される。
［ステップＳ１２４］シミュレーション制御部１２０は、シミュレーション期間内の全タイムステップ分の計算が完了したか否かを判断する。未計算のタイムステップがあれば、処理がステップＳ１０２に進められる。
このようにして、シミュレーション期間内のタイムステップごとの各連続体粒子の位置や物理量が算出される。その結果、解析対象である流体や弾性体の動きが、シミュレートされる。
以上説明したように、第２の実施の形態によれば、連続体粒子が円盤要素の裏側に回り込んだ時に反発力を与えることで、円盤要素の形作る境界形状に細い隙間があっても、連続体粒子の流れ込みを扱うことが可能となる。
また、円盤要素を用いて固定境界の形状を表現しているため、要素間の接続情報を取り扱わずにすむ。これにより、形状を生成するときに煩雑な処理が削減され、表面を隙間なく覆う固定境界の設定が容易となる。
なお円盤要素を用いて一般の形状を扱うと、円盤の端が表現したい形状からはみ出ることがあり得るが、円盤要素同士の角度を用いて、反発力の与え方に修正を施すことで、円盤要素の端において、不自然な反発力が発生することが抑止されている。
第２の実施の形態に示した粒子法を用いたシミュレーションは、広い分野で利用可能である。例えば、河川や海の水の流体解析に対して行うことで防災の計画立案に役立てたり、鋳造過程を解析したりすることで製品設計に用いることができる。弾性体に対して適用することで、製品設計の際に封止ゲルの形状などを適切に決定することができる。特に、非常に狭い隙間に入り込む流体や弾性体の運動を解析するのに、第２の実施の形態に係るシステムは有用である。例えば、弾性体を押しつぶすような動きをする機械により押しつぶされる弾性体の運動を解析する場合、押しつぶされるときの狭い隙間における弾性体の動きを正確に解析することが重要となる。押しつぶす動作により生じる狭い空間に連続体粒子が入り込めないのでは、正確な解析結果は望めないが、第２の実施の形態では狭い隙間に連続体粒子が入り込めるため、正確な解析が可能である。
以上、実施の形態を例示したが、実施の形態で示した各部の構成は同様の機能を有する他のものに置換することができる。また、他の任意の構成物や工程が付加されてもよい。さらに、前述した実施の形態のうちの任意の２以上の構成（特徴）を組み合わせたものであってもよい。
１ 第１の領域
２ 第２の領域
３ 境界
３ａ〜３ｅ 境界要素
４〜７ 粒子
４ａ，７ａ 力
１０ シミュレーション装置
１１ 記憶部
１１ａ 粒子データ
１１ｂ 境界要素データ
１２ 演算部

Claims (5)

  1. コンピュータに、
    連続体を表す複数の粒子の位置と、第１の領域と第２の領域との境界を表し、前記第１の領域を向く第１面と前記第２の領域を向く第２面とを有する複数の要素の位置とに基づいて、前記複数の粒子それぞれを注目粒子としたとき、第２面側に前記注目粒子がある要素が少なくとも１つ存在し、かつ、前記注目粒子の位置が、該要素のうちの前記注目粒子に最も近い最近接要素に対して前記第１の領域側に凸となるように交差する交差要素の第１面側ではない場合、前記注目粒子が前記第２の領域内にあると判断し、第２面側に前記注目粒子がある要素が存在しないか、該要素が存在しても、前記注目粒子の位置が前記交差要素の第１面側にある場合、前記注目粒子は前記第１の領域内にあると判断し、
    前記第２の領域内にある第２粒子に関して、前記複数の要素で表される前記境界と前記第２粒子との距離に基づいて、前記第１の領域側に向かう反発力を算出し、
    前記第１の領域内にある第１粒子には、前記境界と前記第１粒子との距離に基づく、前記境界から離れる方向への力を作用させず、前記第２の領域内にある前記第２粒子には前記反発力を作用させて、前記複数の粒子の運動を解析する、
    処理を実行させるシミュレーションプログラム。
  2. 前記算出では、前記第２粒子が第２面側にある要素のうちの、前記第２粒子に最も近い要素の第１面側の法線方向の前記反発力を算出することを特徴とする請求項１記載のシミュレーションプログラム。
  3. コンピュータに、
    連続体を表す複数の粒子の位置と、第１の領域と第２の領域との境界を表し、前記第１の領域を向く第１面と前記第２の領域を向く第２面とを有する複数の要素の位置とに基づいて、前記複数の粒子それぞれが前記第１の領域内にあるのか前記第２の領域内にあるのかを判断し、
    前記第２の領域内にある第２粒子に関して、前記第２粒子が、前記第１の領域側に凹となるように交差する第１および第２の要素それぞれの第２面側にある特定第２粒子に該当する場合、前記第１および第２の要素のうち、前記第２粒子との距離が長い方の要素における前記第２粒子との距離に基づいて、該要素の第１面側の法線方向の反発力を計算し、前記第２粒子が前記特定第２粒子に該当しない場合、前記複数の要素で表される前記境界と前記第２粒子との距離に基づいて、前記第１の領域側に向かう反発力を算出し、
    前記第１の領域内にある第１粒子には、前記境界と前記第１粒子との距離に基づく、前記境界から離れる方向への力を作用させず、前記第２の領域内にある前記第２粒子に前記反発力を作用させて、前記複数の粒子の運動を解析する、
    処理を実行させるシミュレーションプログラム。
  4. コンピュータが、
    連続体を表す複数の粒子の位置と、第１の領域と第２の領域との境界を表し、前記第１の領域を向く第１面と前記第２の領域を向く第２面とを有する複数の要素の位置とに基づいて、前記複数の粒子それぞれを注目粒子としたとき、第２面側に前記注目粒子がある要素が少なくとも１つ存在し、かつ、前記注目粒子の位置が、該要素のうちの前記注目粒子に最も近い最近接要素に対して前記第１の領域側に凸となるように交差する交差要素の第１面側ではない場合、前記注目粒子が前記第２の領域内にあると判断し、第２面側に前記注目粒子がある要素が存在しないか、該要素が存在しても、前記注目粒子の位置が前記交差要素の第１面側にある場合、前記注目粒子は前記第１の領域内にあると判断し、
    前記第２の領域内にある第２粒子に関して、前記複数の要素で表される前記境界と前記第２粒子との距離に基づいて、前記第１の領域側に向かう反発力を算出し、
    前記第１の領域内にある第１粒子には、前記境界と前記第１粒子との距離に基づく、前記境界から離れる方向への力を作用させず、前記第２の領域内にある前記第２粒子には前記反発力を作用させて、前記複数の粒子の運動を解析する、
    シミュレーション方法。
  5. 連続体を表す複数の粒子の位置と、第１の領域と第２の領域との境界を表し、前記第１の領域を向く第１面と前記第２の領域を向く第２面とを有する複数の要素の位置とを記憶する記憶部と、
    前記複数の粒子の位置と前記複数の要素の位置とに基づいて、前記複数の粒子それぞれを注目粒子としたとき、第２面側に前記注目粒子がある要素が少なくとも１つ存在し、かつ、前記注目粒子の位置が、該要素のうちの前記注目粒子に最も近い最近接要素に対して前記第１の領域側に凸となるように交差する交差要素の第１面側ではない場合、前記注目粒子が前記第２の領域内にあると判断し、第２面側に前記注目粒子がある要素が存在しないか、該要素が存在しても、前記注目粒子の位置が前記交差要素の第１面側にある場合、前記注目粒子は前記第１の領域内にあると判断し、前記第２の領域内にある第２粒子に関して、前記複数の要素で表される前記境界と前記第２粒子との距離に基づいて、前記第１の領域側に向かう反発力を算出し、前記第１の領域内にある第１粒子には、前記境界と前記第１粒子との距離に基づく、前記境界から離れる方向への力を作用させず、前記第２の領域内にある前記第２粒子には前記反発力を作用させて、前記複数の粒子の運動を解析する演算部と、
    を有するシミュレーション装置。
JP2015000797A 2015-01-06 2015-01-06 シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置 Active JP6458501B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2015000797A JP6458501B2 (ja) 2015-01-06 2015-01-06 シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置
US14/950,373 US10068037B2 (en) 2015-01-06 2015-11-24 Simulation method and simulation apparatus for continuum motion analysis using a particle method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2015000797A JP6458501B2 (ja) 2015-01-06 2015-01-06 シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置

Publications (2)

Publication Number Publication Date
JP2016125934A JP2016125934A (ja) 2016-07-11
JP6458501B2 true JP6458501B2 (ja) 2019-01-30

Family

ID=56286668

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2015000797A Active JP6458501B2 (ja) 2015-01-06 2015-01-06 シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置

Country Status (2)

Country Link
US (1) US10068037B2 (ja)
JP (1) JP6458501B2 (ja)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6547547B2 (ja) * 2015-09-25 2019-07-24 富士通株式会社 粒子シミュレーションプログラム、計算機資源配分方法、および粒子シミュレーション装置
JP6458776B2 (ja) 2016-06-24 2019-01-30 株式会社デンソー 車両用表示装置
JP6697407B2 (ja) * 2017-03-07 2020-05-20 公益財団法人鉄道総合技術研究所 流体シミュレーション方法及び流体シミュレーションのプログラム
JP6897477B2 (ja) * 2017-10-10 2021-06-30 富士通株式会社 流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置
CN111104728B (zh) * 2019-10-15 2023-04-25 江汉大学 结构表面磨损的仿真预测方法和装置
JP7285223B2 (ja) * 2020-01-21 2023-06-01 住友重機械工業株式会社 シミュレーション装置、シミュレーション方法、及びプログラム

Family Cites Families (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140226158A1 (en) * 2004-03-06 2014-08-14 Michael Trainer Methods and apparatus for determining particle characteristics
JP4798661B2 (ja) 2006-10-27 2011-10-19 みずほ情報総研株式会社 流体解析装置、流体解析方法及び流体解析プログラム
KR100984048B1 (ko) * 2008-10-14 2010-09-30 한국전자통신연구원 파티클 유체 시뮬레이션에서의 강성체 상호작용 처리 방법
KR101159163B1 (ko) * 2008-12-10 2012-06-22 한국전자통신연구원 비압축 상태의 유체 시뮬레이션 방법, 기록매체, 장치
JP2010243293A (ja) * 2009-04-03 2010-10-28 Toshiba Corp 流動解析方法、流動解析装置、及び流動解析プログラム
US8762117B2 (en) * 2009-06-03 2014-06-24 Algoryx Simulation Ab Method, an apparatus and computer program product for simulating dynamic fluids
WO2012025995A1 (ja) * 2010-08-24 2012-03-01 富士通株式会社 連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置
JP5670832B2 (ja) * 2011-05-24 2015-02-18 富士通株式会社 シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
US8903693B2 (en) * 2011-06-24 2014-12-02 Siemens Aktiengesellschaft Boundary handling for particle-based simulation
JP5844165B2 (ja) * 2012-01-26 2016-01-13 住友重機械工業株式会社 解析装置
JP6048062B2 (ja) * 2012-10-18 2016-12-21 富士通株式会社 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
JP2014104696A (ja) * 2012-11-29 2014-06-09 Japan Steel Works Ltd:The 演算プログラム及び演算方法
US20150161810A1 (en) * 2013-12-10 2015-06-11 Nvidia Corporation Position based fluid dynamics simulation

Also Published As

Publication number Publication date
JP2016125934A (ja) 2016-07-11
US20160196373A1 (en) 2016-07-07
US10068037B2 (en) 2018-09-04

Similar Documents

Publication Publication Date Title
JP6458501B2 (ja) シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置
Shojaei et al. An adaptive multi-grid peridynamic method for dynamic fracture analysis
Idelsohn et al. A Lagrangian meshless finite element method applied to fluid–structure interaction problems
Ingram et al. Developments in Cartesian cut cell methods
Ferrari et al. A new 3D parallel SPH scheme for free surface flows
Janßen et al. On enhanced non-linear free surface flow simulations with a hybrid LBM–VOF model
Edwards et al. Detailed water with coarse grids: Combining surface meshes and adaptive discontinuous galerkin
JP2009069930A (ja) 粒子法シミュレーションのためのスライスデータ構造、およびスライスデータ構造を利用した粒子法シミュレーションのｇｐｕへの実装方法
US8040347B2 (en) Method for constructing surface of fluid-body simulation based on particle method, program for the same, and storage medium for storing program
US10503561B2 (en) Particle simulation apparatus and computer resource allocating method
Samulyak et al. Lagrangian particle method for compressible fluid dynamics
JP5704246B2 (ja) 物体運動解析装置、物体運動解析方法、及び物体運動解析プログラム
US20200110911A1 (en) Particle simulation device, particle simulation method, and particle simulation program
Quan et al. Meshing molecular surfaces based on analytical implicit representation
Samuel A deformable particle-in-cell method for advective transport in geodynamic modelling
Bähr et al. Stable honeycomb structures and temperature based trajectory optimization for wire-arc additive manufacturing
KR100588000B1 (ko) 유체 애니메이션에서의 자유경계 추적 장치 및 그 방법
Reuther et al. A comparison of methods for the calculation of interface curvature in two-dimensional cellular automata solidification models
Staroszczyk Simulation of solitary wave mechanics by a corrected smoothed particle hydrodymamics method
JP5454693B2 (ja) 連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置
KR102520732B1 (ko) 유동해석 데이터 처리장치 및 그 장치에서 각 기능을 실행시키기 위해 매체에 저장된 컴퓨터 프로그램
Sandim et al. Boundary particle resampling for surface reconstruction in liquid animation
JP2022041425A (ja) シミュレーションプログラム、シミュレーション方法およびシミュレーションシステム
Lu et al. Parallel contact-aware simulations of deformable particles in 3D Stokes flow
US11835054B2 (en) Method for automatic detection of axial cooling fan rotation direction

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20171113

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20180427

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20180522

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20180702

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20181127

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20181210

R150 Certificate of patent or registration of utility model

Ref document number: 6458501

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150