JP5454693B2 - 連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置 - Google Patents

Description

本発明は、連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置に関する。

数値計算を用いて流体である水や空気の流れを解析する流体解析や、圧縮された弾性体であるゴムの振る舞いを解析する弾性体解析などといった連続体の運動の解析において、粒子法と呼ばれる技術が提案されている。具体的には、粒子法とは、連続体の運動を有限の数の粒子の運動として解析する手法である。現在提案されている代表的な粒子法としては、SPH(Smoothed Particles Hydrodynamics)法やMPS(Moving Particles Semi-implicit)法といったものがある。また、粒子法には、粉流体を粒子の集まりとして表現し粉流体の挙動を解析する手法である個別要素法と呼ばれる方法なども含まれる。

従来、粒子法において、粒子分布を用いて解析対象を表現する場合の固定境界を表現する方法としては、以下のような方法が提案されてきた。

まず、粒子法の標準的な技術として以下の三つの技術が提案されている。まず第１の方法としては、境界に沿って仮想的な流体粒子（以下、「境界粒子」という。）を配置する。そして、境界粒子を固定境界に固着した流体とみなし、連続体を表現する粒子と境界粒子の間に働く力を一般的な流体解析において働く力として与えることにより、固定境界を表現する技術が提案されている。第２の方法は、第１の方法と同様に境界粒子を配置する。そして、各粒子の間の距離に基づく関数で大きさが与えられ、且つ各粒子の相対位置ベクトルに沿った方向に働く、固定境界から各粒子への反発力を設定して、連続体を表現する技術が提案されている。この技術では、連続体を表現する粒子運動の領域を固定境界の内側の領域に限定することにより解析を行うことができる。第３の方法は、解析対象である連続体を表現する粒子に対して、予め設定された境界面に対して鏡映対称な位置に仮想的な連続体粒子を配置することで固定境界を表現する技術が提案されている。

また、他の方法としては、まず境界粒子及び連続体粒子のそれぞれが球面とは限らない滑らかな閉曲面を有するものとみなす。そして、それら閉曲面が重なりを持ったときに適宜縮小し、閉曲面が一点で接するときに、各粒子中心から当該接点までの距離と方向に基づいて固定境界から各粒子に働く反発力などを計算することで固定境界を表現する技術が提案されている。

また、他の方法としては、まず境界粒子を配置する。そして、境界粒子との距離が一定以下に近づいた場合の連続体粒子において、連続体粒子から下ろした垂線の距離及び境界面に沿った境界粒子と連続体粒子との距離を求める。そして、求めたそれぞれの距離による関数で大きさが与えられ、且つ境界面の法線に沿った向きを有する、固定境界から各粒子への反発力を、連続体粒子に対して与えて、固定境界を表現する技術が提案されている。

さらに、他の方法としては、まず固定境界をポリゴンの集合として表現する。そして、連続体粒子は重心座標が最も近いポリゴンとの間に、連続体粒子とポリゴンを含む平面との距離の関数で表される反発力を与える技術が提案されている。この技術では、連続体粒子とポリゴンを含む平面との距離が一定値ε以下になる場合にのみ、平面から連続体粒子に反発力を及ぼすものとしている。

特開２００９−０２６２７９号公報

G.R. Liu, M.B. Liu 著 「Smoothed Particle Hydrodynamics: A Meshfree Particle Method」World Scientific Pub Co Inc 出版 2003年"p.138-141 4.4.8節 Boundary treatment" M.G. Gesteira, B.D. Rogers, R.A. Dalrymplem, A.J.C. Crespo, M. Narayanaswamy 著 「User Guide for the SPHysics Code v1.4」"p.16-19 2.4節 Boundary conditions（ただし、2.4.3節を除く）" 原田隆宏、越塚誠一、河口洋一郎 著 「流体と布のリアルタイム連成シミュレーション」 情報処理学会研究報告、グラフィクスとＣＡＤ研究会報告 2007/11/12（全文）

しかしながら、一般的な流体解析において境界粒子から働く力により固定境界を表現する従来技術では、滑らかな形状の固定境界を表現するためには多数の粒子を用いる必要があり、解析を行う計算機のメモリの使用量及び計算時間のいずれも多く必要としていた。これは、固定境界からの反発力を粒子間の相対位置ベクトルに沿った向きに設定する従来技術、及び境界面に対して鏡映対称な位置に仮想的な連続体粒子を配置する技術のいずれにおいても同様である。

加えて、境界面に対して鏡映対称な位置に仮想的な連続体粒子を配置する技術では、鏡映対称な位置を計算するために境界面の解析的に表現しなくてはならない。そのため、法線方向を表す数式が特に与えられていない又は与えるのが困難な境界形状に対しては適用できないという問題があった。

また、任意の滑らかな閉曲面を連続体粒子、境界粒子の両方に設定する従来技術では、固定境界からの反発力を算出するためには連続体粒子及び境界粒子の対毎に多くのデータ及び計算量が必要となる。そのため、この技術でもメモリの使用量及び計算時間を多く必要としていた。

これに対して、固定境界からの力として法線方向の反発力を設定する従来技術では、メモリの使用量などを減らすことはできる。しかし、この従来技術は、ある境界粒子ｉからｘ座標の方向に正又は負の方向に一つずれ、さらにｙ座標の方向に正又は負の方向に一つずれた位置にある４つの近傍粒子（ｉ−１、ｉ＋１、ｊ−１、ｊ＋１）を求める。そして、ｉ−１からｉ＋１へ向かうベクトルとｊ−１からｊ＋１へ向かうベクトルの外積により法線ベクトルを求める。すなわち、この従来技術では、境界面上の粒子分布が正則な格子点又はそれらを連続的に変形させた座標位置に配置する必要がある。このため、例えばテトラメッシュの重心位置に境界粒子を配置するなどした場合には応用ができず、適用できない場合があるという問題があった。

また、重心座標が最も近いポリゴンとの間に反発力を与える従来技術では、ポリゴンのサイズが一定でなく異なるポリゴンのサイズを用いた場合には、ポリゴンの境界付近においてポリゴンのサイズの不均一に寄る不合理な反発力を生じるおそれがあり、適用するにはポリゴンのサイズが一定である必要がある。そのため、この従来技術も適用できない場合があるという問題があった。

本発明は、上記に鑑みてなされたものであって、適用できない場合を減少させるとともに、解析を行うために必要なメモリ容量及び計算時間を軽減する連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置を提供することを目的とする。

本願の開示する連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置は、一つの態様において、連続体を粒子として表した場合の各連続体粒子を表す影響範囲の情報を含むデータ、及び固定境界を任意の形状を有する微小な領域の集合として表した場合の各微小面要素を表す法線方向の情報を含むデータの入力を受け付ける手順と、各前記連続体粒子の前記影響範囲に含まれる前記微小面要素から各前記連続体粒子に対して働く前記微小面要素の前記法線方向の反発力を算出する手順と、各前記反発力を足し合わせることで各前記連続体粒子に対して働く前記固定境界からの力を求める手順と、各前記連続体粒子に対して働く前記固定境界からの力を基に、前記連続体の運動の解析を行う手順とをコンピュータに実行させる。

本願の開示する連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置の一つの態様によれば、適用できない場合を減少させるとともに、解析を行うために必要なメモリ容量及び計算時間を軽減するという効果を奏する。

図１は、本実施例に係る連続体運動解析装置のブロック図である。 図２は、本実施例に係る連続体運動解析プログラムによる連続体の運動の解析処理のフローチャートである。 図３は、連続体粒子と微小面要素との相対位置ベクトルの微小面要素の法線方向の成分と法線に直交する方向の成分とを説明するための図である。 図４は、連続体粒子の影響範囲と微小面要素とが交わるか否かの判定を説明するための図である。 図５は、連続体粒子と微小面要素との共通部分の面積を説明するための図である。 図６は、連続体粒子に接している微小面要素のリストを説明するための図である。 図７は、連続体粒子方向に突出するように交わる微小面要素の図である。 図８は、最近接微小面要素を同定する処理のフローチャートである。 図９は、２つの微小面要素が連続体粒子方向に突出するように交わるか否かの判定を説明するための図である。 図１０は、２つの微小面要素が交わるか否かを判定する処理のフローチャートである。 図１１は、一方の微小面要素を含む平面と他方の微小面要素とが交わっている状態の図である。 図１２は、連続体運動解析プログラムを実行するコンピュータを示す図である。

以下に、本実施例に係る連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置の実施例を図面に基づいて詳細に説明する。なお、以下の実施例により本願の開示する連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置が限定されるものではない。

ここでは、コップへの水の注入やダムからの水の流れというような固定境界に対する連続体の運動を解析する場合について説明する。また、以下では、解析対象である連続体を粒子として表したときの粒子を「連続体粒子」と言う。また、固定境界を微小面の集合で表したときの微小面のことを「微小面要素」と言う。連続体粒子とは、言い換えれば、連続体を粒子によってモデル化したものである。また、微小面要素とは、言い換えれば、流体を格納する容器の壁といった固定境界を微小部分に分割した面要素をモデル化したものである。微小面要素は、任意の形状を有する。すなわち、微小面要素は、特定の形状に固定されず自由な形状を持つことができる。具体的には、微小面要素は、固定境界を覆うように配置された微小な円盤で表されたり、同じく固定境界を覆うように配置されたポリゴンとして表わされたりする。さらに、連続体の運動を解析するには予め決められた所定時間の中で連続体の運動を求める必要があり、そのためには、所定時間内の所定の時刻における連続体の運動を求めそれらを一連の運動として統合することになる。

図１は、本実施例に係る連続体運動解析装置のブロック図である。図１に示すように、本実施例に係る連続体運動解析装置は、ユーザインタフェース１、記憶部２、反発力算出部３、圧力算出部４、解析部５及び出力部６を有している。

ユーザインタフェース１は、キーボードやマウスなどの入力装置及びモニタなどの表示装置を備えている。操作者は、ユーザインタフェース１を用いて連続体の運動の解析に用いるデータの入力を行うとともに、表示装置に表示された解析結果を確認する。

ここで、連続体の運動の解析に用いられるデータとは、連続体粒子に設定された物理モデル、各連続体粒子を表現する連続体粒子の情報及び各微小面要素を表現する微小面要素の情報である。ここで、物理モデルとは、連続体の運動の物理表現を表現するためのモデルである。例えば、連続体粒子に設定する物理モデルとしては、流体をモデル化した場合の物理モデルなどがある。そして、連続体粒子の情報とは、例えば、中心座標、速度、影響半径、密度、質量、変形勾配テンソル、材料物性及び温度などである。ここで、連続粒子の中心を中心とし影響半径を半径とする球は、連続体粒子が圧力を加えるなどの影響を他の粒子に与え、また他の粒子から連続体粒子に対して圧力を加えられるなどの影響が与えられる範囲を表す。すなわち、この連続体粒子の中心を中心とし影響半径を半径とする球の外にある粒子はその球の中心にある連続体粒子に対して影響を与えない。また、微小面要素を円盤とした場合、微小面要素の情報とは、例えば、中心座標、微小面要素の法線方向を表す法線ベクトル及び面積などである。ここで、法線方向とは、固定境界から連続体が存在する側に向いた微小面要素における法線の方向のことを指す。また、微小面要素をポリゴンとした場合、微小面要素の情報とは、例えば、各微小面要素の頂点の座標である。ここで、微小面要素をポリゴンとした場合、各微小面要素の頂点の座標が与えられれば、その面積及び法線方向が算出できる。すなわち、微小面要素をポリゴンとした場合、微小面要素の各頂点の座標を与えることは、微小面要素の面積及び法線方向を与えることと同じである。

記憶部２は、ハードディスクなどの記憶装置である。記憶部２は、ユーザインタフェース１から入力された連続体粒子の情報及び微小面要素の情報を記憶する。

反発力算出部３は、連続体粒子の影響範囲内にある微小面要素を抽出する。そして、反発力算出部３は、抽出した各微小面要素からの微小面要素の法線方向に向けた連続体粒子にかかる力を算出し、さらに連続体粒子と微小面要素の影響範囲との交わる部分の面積に応じた重み付けを行って、各微小面要素による反発力を算出する。この反発力の算出は後で詳細に説明する。

圧力算出部４は、連続体粒子毎に反発力算出部３が算出した反発力を加算し、各連続体粒子に対する固定境界面からの圧力を算出する。

解析部５は、ある時刻における各連続体粒子における物理モデルに基づいた加速度を算出する。また、解析部５は、各連続体粒子にかかる固定境界からの圧力による加速度を計算する。そして、解析部５は、連続体粒子の物理モデルに基づいた加速度に、固定境界からの圧力による加速度を、加算することにより連続体粒子の加速度を求める。

そして、解析部５は、求めた加速度を含む時間微分項を用いて連続体粒子の物理量に対して時間積分し、ある時刻における連続体粒子の運動を取得する。解析部５は、このある時刻における連続体粒子の運動を連続する各時刻について行い、連続体の運動の解析を行う。

出力部６は、解析部５による解析結果を出力する。例えば、出力部６は、解析結果を記憶装置にファイルとして出力したり、表示装置に図として表示したり、又はプリンタ装置を用いて紙に印刷して出力したりする。

以上の説明のように、まず、反発力算出部３により、連続体粒子の影響範囲に近づいた微小面要素の法線方向からの反発力が求められる。そして、圧力算出部４により、連続体粒子の影響範囲に近づいた微小面要素からの連続体粒子に対する圧力が求められる。そして、解析部５により、連続体粒子に対する圧力から加速度が求められる。これにより、微小面要素が影響範囲以内に近づいた連続体粒子には、その微小面要素の法線方向の加速度が加えられることになり、連続体粒子に対する固定境界からの作用を表す固定境界条件を連続体の運動に適用することができる。次に、固定境界条件を連続体に適用して連続体の運動を解析する具体的な処理について詳細に説明する。

まず、図２を参照して、本実施形態に係る連続体運動解析プログラムによる連続体の運動の解析処理の流れについて説明する。図２は、本実施例に係る連続体運動解析プログラムによる連続体の運動の解析処理のフローチャートである。

記憶部２は、ユーザインタフェース１から入力された連続体粒子の情報及び微小面要素の情報を取得し記憶する（ステップＳ１）。

解析部５は、入力された連続粒子に設定された物理モデルに従って連続粒子の運動を解き、物理モデルに基づく各連続体粒子の加速度を算出する（ステップＳ２）。ここで、解析部５は、運動解析を行うための一般的な運動方程式を用いて解析を行う。例えば、解析部５は、流体をモデル化した連続体粒子に対しては、ナヴィエ・ストークス方程式を解いて、加速度や温度の時間微分を求める。ナヴィエ・ストークス方程式は非圧縮性の流れの場合には次の数式１のように表される方程式である。

ここで、ｖは流速、Ｆは単位質量あたりに働く外力、ρは密度、ｐは圧力、νは動粘性計数である。

具体的には、解析部５は、ナヴィエ・ストークス方程式及び質量保存則から導かれる連続の式、さらに必要であれば状態方程式などを連立して数値的に解くことで解析を行う。

次に、反発力算出部３は、連続体粒子のそれぞれに対して、影響範囲内にある微小面要素のリストを作成する（ステップＳ３）。ここで、図３〜図６を参照して、反発力算出部３による、連続体粒子の影響範囲内にある微小面要素のリスト作成について詳細に説明する。ここで、図３は、連続体粒子と微小面要素との相対位置ベクトルの微小面要素の法線方向の成分と法線に直交する方向の成分とを説明するための図である。また、図４は、連続体粒子の影響範囲と微小面要素とが交わるか否かの判定を説明するための図である。また、図５は、連続体粒子と微小面要素との共通部分の面積を説明するための図である。また、図６は、連続体粒子に接している微小面要素のリストを説明するための図である。

以下の説明では、連続体粒子すべてに番号をふり、任意の連続体粒子の番号をiとする。また、微小面要素すべてに番号をふり、任意の微小面要素の番号をjとする。そして、番号iの連続体粒子を連続体粒子iと表し、番号jの微小面要素を微小面要素jと表す。さらに、図３に示すように、連続体粒子iと微小面要素jとの相対位置ベクトルをr_ijとする。ここで、連続体粒子iと微小面要素jとの相対位置ベクトルとは、連続体粒子iの中心から微小面要素jの中心へと向かうベクトルである。また、連続体粒子iの影響半径をh_iとする。さらに、微小面要素jの面積をA_jとする。図３では、説明が分かり易いように、連続体粒子の影響範囲にiという符号をふっているが、iは連続体粒子を表しているものとする。

反発力算出部３は、図４に示すr_ijのn_j方向成分r_nijを、r_nij=(r_ij・n_j) n_jとして求める。また、反発力算出部３は、図４に示すr_ijのn_jに垂直な成分r_tijを、r_tij=r_ij-r_nijとして求める。また、反発力算出部３は、図４に示す微小面要素の半径R_jを、R_j=(A_j/π)^1/2として求める。

ここで、|r_nij|＜h_iのとき、微小面要素jを含む面と連続体粒子iの影響範囲とが交わっている場合には、円盤Ｄを構成する。この円盤Ｄの半径R_pijは、R_pij=(h_i ²-r_nij ²)^1/2として表される。

そして、図４に示すように、連続体粒子iの影響範囲と微小面要素jが交わっている場合には、R_jとR_pijとを加えた長さよりもr_tijの長さのほうが短くなる。すなわち、反発力算出部３は、|r_nij|＜h_i 且つR_j+R_pij＞|r_tij|という条件が満たされるときに、微小面要素jは連続体粒子iの影響範囲と交わっている、言い換えれば、微小面要素jは連続体粒子iの影響範囲内にあると判定する。

そして、反発力算出部３は、微小面要素jが連続体粒子iの影響範囲と交わっているか否かの判定を、全ての微小面要素j及び連続体粒子iの組について行うことで、各連続体粒子iの影響範囲と交わっている微小面要素jを抽出する。

さらに、反発力算出部３は、図５に示す連続体粒子iの影響範囲と微小面要素jが交わっている場合の交わっている部分の面積S_ijを求める。具体的には、反発力算出部３は、中心間の距離が|r_tij|離れている半径R_jの円と半径R_pijの円の共通部分の面積として面積S_ijを計算する。

そして、反発力算出部３は、連続体粒子i毎に、連続体粒子iの影響範囲と交わる微小面要素jを表す番号、相対位置ベクトルの法線成分及び法線に垂直な成分、並びに連続体粒子iの影響範囲と微小面要素jが交わる部分の面積を記載したリストを作成する。このリストを概念的に表したものが図６である。図６に示すように、各連続体粒子に対してリストが作成される。例えば、図６では、連続体粒子i-1、i、i+1、i+2に対するリストを表示している。リストには、枠で囲われたグループ１００が一つの微小面要素のグループとして記載される。そして、グループ１００には、欄１０１に微小面要素を表す番号、欄１０２に相対位置ベクトルの法線成分、欄１０３には相対位置ベクトルの法線に垂直な成分、及び欄１０４には連続体粒子の影響範囲と微小面要素が交わる部分の面積が記載される。例えば、連続体粒子iの影響範囲に連続体要素jが交わっている場合には、グループ１００の欄１０１には微小面要素番号j、欄１０２にはr_tij、欄１０３にはr_nij、欄１０４にはS_ijが記載される。また、図６の連続体粒子iに対するリストのように、微小面要素番号k、微小面要素番号jというように、微小面要素毎に微小面要素に対応するグループ１００が作成される。また、図６では、連続体粒子i-1、i、i+1、i+2のいずれに対するリストにも微小面要素番号jが含まれているので、連続体粒子i-1、i、i+1、i+2のいずれに対しても微小面要素番号jが影響範囲内にあることが分かる。

リスト完成後、反発力算出部３は、リストに存在する微小面要素jと連続体粒子iの影響範囲との交わり部分の面積の総和を求める（ステップＳ４）。ここで、微小面要素と連続体粒子の影響範囲との交わり部分の面積の総和は、例えば、連続体粒子が連続体粒子iの場合であれば、総和S_total,iと表すことにする。

次に、反発力算出部３は、各微小面要素からそれぞれの連続体粒子に働く反発力を算出する（ステップＳ５）。以下では、反発力算出部３による反発力の算出について詳細に説明する。ここでは、微小面要素jから連続体粒子iに対して働く反発力を例に説明する。

反発力算出部３は、微小面要素jの法線に沿った連続体粒子iと微小面要素jとの距離|r_nij|が、一定距離r_cutoffに近づくと反発力が働き、距離|r_nij|が短くなると反発力は大きくなり、|r_nij|が０になる極限でも発散しない反発力を算出する式を予め記憶している。例えば、本実施例では、反発力算出部３は次のような数式２を記憶している。

ここで、r_fmaxは、反発力が距離に反比例して大きくなる値をとることのできる距離の限界値として決定した値であり、０に近い固定値である。すなわち、連続体粒子ｉと微小面要素ｊとの距離がr_fmaxよりも近づくと反発力は一定になる。r_fmaxは、０に近い値であれば任意の値を取ることができる。ただし、rfmaxを大きくしすぎると十分な反発力が得られなくなり、逆にr_fmaxを小さくしすぎると反発力が大きくなりすぎて数値計算が安定しなくなる。そこで、r_fmaxは、連続体を解析するに当たり要求される反発力と数値計算の安定を考慮して設定されることが好ましい。統計的には、r_fmaxをr_cutoffの１００分の１より小さくした場合、数値計算が安定しなくなる。また、統計的には、r_fmaxをr_cutoffの１０分の１より大きくした場合、十分な反発力が得られなくなり、固定境界が不明になってしまう。そこで、本実施例ではr_fmaxは、r_cutoffの１０分の１〜１００分の１としている。このr_fmaxが所定距離の一例にあたる。また、q=|r_nij|/r_cutoff,q_fmax=r_fmax/r_cutoffである。すなわち、qはr_cutoffで規格化された距離を表す値である。そして、q_fmaxは、反発力が距離に反比例して大きくなる値をとることのできる距離の限界値であるr_fmaxをr_cutoffで規格化した値であり、固定値となる。また、f_0,iは基準となる反発力の大きさである。f_0,iは大きな値であれば特に制限は無い。例えば、本実施例では、f_0,i=m_i×c_x,i ²/h_i（m_iは連続体粒子iの質量、c_x,iは音速）としている。

数式２では、|r_nij|＜r_fmaxの場合、最後の項が(1-q_fmax)/q_fmaxとなっている。ここで、q_fmaxは固定値であるので、(1-q_fmax)/q_fmaxも固定値となる。すなわち、|r_nij|＜r_fmaxの場合、r_nijがいくら短くなっても値は増えず、r_nijが０になっても反発力の値は発散しない。ここで、r_nijが０になっても発散しないように数式２が構成されている理由について、図７を参照して説明する。図７は、連続体粒子方向に突出するように交わる微小面要素の図である。

図７の斜線で表される部分が固定境界に対して連続体と反対側の領域である。そして、図７では微小面要素２０１と微小面要素２０２とが形成する面が連続体粒子方向に突出するように、微小面要素２０１と微小面要素２０２とが交わっている。以下では、１つの連続体粒子２１０が矢印２１１に沿って動く場合で説明する。

微小面要素２０１及び微小面要素２０２から連続体粒子２１０に対して働く反発力は、斜線で表される領域とは反対に向く法線方向に働く。矢印２１１のように連続体粒子２１０が移動する場合、微小面要素２０２と矢印２１１との交点２１２を越えるまでは、連続体粒子２１０には微小面要素２０２からの反発力は働かず、微小面要素２０１からの反発力のみが働く。そして、連続体粒子２１０が交点２１２を通過する瞬間に、微小面要素２０２と連続体粒子２１０の距離は０となる。ここで、もし連続体粒子と微小面要素との距離が０の場合に発散する計算式を用いた場合、連続体粒子２１０が交点２１２を通過する瞬間に、極端に大きな力が微小面要素２０２から連続体粒子２１０に対して働いてしまう。そのため、連続体粒子と微小面要素との距離が０の場合に発散する式を用いた場合、図７のような状況下では反発力の計算が破綻してしまう。そこで、図７のような場合にも計算が破綻しないように、数式２はr_nijが０になっても発散しないように構成されている。こうすることにより、数値計算を安定的に進めることができる。

そして、反発力算出部３は、リストの全ての微小面要素について反発力の算出が完了したか否かを判定する（ステップＳ６）。リストの全ての微小面要素について反発力の算出が完了していない場合（ステップＳ６否定）、ステップＳ５に戻る。

これに対して、リストの全ての微小面要素について反発力の算出が完了している場合（ステップＳ６肯定）、反発力算出部３は、各連続体粒子に対する最近接微小面要素の同定を行う（ステップＳ７）。この最近接微小面要素を同定する処理については後で詳述する。

反発力算出部３は、連続体粒子iに対する最近接微小面要素からの反発力を、距離０で発散するよう定式化された反発力f_nearest,ijに変更する（ステップＳ８）。本実施例では、数式３で表される反発力を最近接微小面要素である微小面要素jからの連続体粒子iに対する反発力として設定する。

このように、最近接微小面要素からの反発力を距離０で発散するように設定することで、きわめて大きな運動エネルギーを持った連続体粒子であっても固定境界を突き抜けるという事態を回避することができる。

圧力算出部４は、数式２で表される連続体粒子iへの反発力f_ij及び数式３で表される連続体粒子iへの反発力f_nearest,ijを加算して、連続体粒子iへの固定境界からの圧力を算出する（ステップＳ９）。

解析部５は、固定境界からの圧力による連続体粒子iの加速度を算出する。そして、解析部５は、固定境界からの圧力による連続体粒子iの加速度と物理モデルを基に求めた連続体粒子iの加速度とを加算する（ステップＳ１０）。

解析部５は、全ての連続体粒子についてその時刻における加速度の計算が完了したか否かを判定する（ステップＳ１１）。加速度を計算していない連続体粒子がある場合（ステップＳ１１否定）、ステップＳ３に戻る。

そして、全ての連続体粒子についてその時刻における加速度の計算を完了した場合（ステップＳ１１肯定）、解析部５は、各連続体粒子の加速度を含む時間微分項を用いて各連続体粒子の物理量を時間積分し、その時刻における各連続体粒子の運動を求める（ステップＳ１２）。ここで、「各連続体粒子の運動を求める」とは、各連続体粒子の中心座標、速度、影響半径、密度、質量、変形勾配テンソル、材料物性及び温度等の物理量を求めることを言う。

解析部５は、運動を求めた時刻における連続体粒子の運動の計算結果を出力する（ステップＳ１３）。

解析部５は、予め決められた所定時間内における全てのタイミングでの連続体粒子の運動の計算が完了したか否かを判定する（ステップＳ１４）。所定時間内における全てのタイミングでの連続体粒子の運動の計算が完了していない場合（ステップＳ１４否定）、ステップＳ２に戻る。その後、反発力算出部３、圧力算出部４及び解析部５は、時間を進めて次のタイミングにおける連続体粒子の運動の計算を行う。

これに対して、所定時間内における全てのタイミングでの連続体粒子の運動の計算が完了した場合（ステップＳ１４肯定）、解析部５は、各時刻での連続体粒子の運動をまとめて、所定時間内における連続体粒子全体の運動を求める（ステップＳ１５）。これにより、所定時間内における連続体の一連の運動が求まり、連続体運動解析装置は、連続体の運動の解析を終了する。

次に、最近接微小面要素を同定する処理について、図８を参照して詳細に説明する。ここでの説明は、図２のステップＳ７における、反発力算出部３による各連続体粒子に対する最近接微小面要素を同定する処理の詳細な説明である。図８は、最近接微小面要素を同定する処理のフローチャートである。以下では、連続体粒子iに対する最近接微小面要素の同定について説明する。

反発力算出部３は、作成した連続体粒子iの影響半径と交わる微小面要素のリストの中から、次の条件を満たす微小面要素jを抽出する（ステップＳ７１）。この条件とは、連続体粒子iが微小面要素jの法線方向に位置し、且つ連続体粒子iから微小面要素jに対して垂線を下ろすと、微小面要素j内に接点を持つ、という２つの条件である。具体的には、反発力算出部３は、連続体粒子iに対して微小面要素jとがr_ij・n_j＜０の関係を満たすか否かを判定する。「・」の記号は内積を表すものとする。ここで、図３で示されるように、r_ij・n_j＜０であれば、連続体粒子iが微小面要素jの法線方向に位置しているといえる。さらに、反発力算出部３は、連続体粒子iと微小面要素jとがR_j＞|r_tij|を満たすか否かを判定する。ここで、図４で示されるように、R_j＞|r_tij|であれば、連続体粒子iから微小面要素jに対して垂線を下ろすと、微小面要素j内に接点を持つといえる。そして、反発力算出部３は、r_ij・n_j＜０かつR_j＞|r_tij|を満たす微小面要素jを抽出する。

そして、反発力算出部３は、抽出した微小面要素jの中から連続体粒子iに最も近い微小面要素、すなわち|r_nij|が最も小さくなる微小面要素jを仮の最近接微小面要素とする（ステップＳ７２）。

次に、反発力算出部３は、仮の最近接微小面要素に対して連続体粒子方向に突出するように交わる微小面要素があるか否かを判定する（ステップＳ７３）。連続体粒子方向に突出するように交わる微小面要素があるか否かを判定する方法は後で詳細に説明する。

そして、連続体粒子方向に突出するように交わる微小面要素がある場合（ステップＳ７３肯定）、反発力算出部３は、最近接微小面要素を特定する対象となっている微小面要素の中から、仮の最近接微小面要素として設定した微小面要素を除く（ステップＳ７４）。そして、ステップＳ７２に戻り、反発力算出部３は、再度連続体粒子に最も近い微小面要素を探す。

これに対して、連続体粒子方向に突出するように交わる微小面要素が無い場合（ステップＳ７３否定）、反発力算出部３は、仮の最近接微小面要素としている微小面要素を最近接微小面要素として同定する（ステップＳ７５）。

以上のようにして、反発力算出部３は、図２のステップＳ７における、各連続体粒子に対する最近接微小面要素の同定を行うことができる。

ここで、２つの微小面要素がある場合に、それら２つが連続体粒子方向に突出するように交わるか否かを判定する方法について、図９、図１０及び図１１を参照して説明する。この説明は、図８のステップＳ７３の処理の詳細な説明である。図９は、２つの微小面要素が連続体粒子方向に突出するように交わるか否かの判定を説明するための図である。また、図１０は、２つの微小面要素が交わるか否かを判定する処理のフローチャートである。また、図１１は、一方の微小面要素を含む平面と他方の微小面要素とが交わっている状態の図である。

まず、図９を参照して、２つの微小面要素が連続体粒子方向に突出するように交わっているか否かの判定について説明する。ここでは、２つの微小面要素の番号をj1及びj2として説明する。２つの微小面要素は、微小面要素j1及び微小面要素j2として表される。そして、微小面要素j1の法線ベクトルをn_j1、半径をR_j1とする。また、微小面要素j2の法線ベクトルをn_j2、半径をR_j2とする。さらに、微小面要素j1の中心から微小面要素j2の中心へと向かう相対位置ベクトルをr_j1j2とする。

そして、r_j1j2・n_j1＜０及び-r_j1j2・n_j2＜０の両方の式を満たすとき、微小面要素j1と微小面要素j2とは連続体粒子方向に突出しているといえる。したがって、反発力算出部３は、微小面要素j1を仮の最近接微小面要素としている場合、r_j1j2・n_j1＜０及び-r_j1j2・n_j2＜０の両方の式を満たす微小面要素j2を抽出する。

次に、図１０及び図１１を参照して、反発力算出部３による、微小面要素j1と微小面要素j2とが交わっているか否かの判定方法を説明する。

まず、反発力算出部３は、微小面要素j1を含む平面P1と微小面要素j2が交わっているか否かを判定する（ステップＳ７３１）。具体的には、反発力算出部３は、|r_j1j2・n_j1|＜R_j2sinθであるか否かを判定し、|r_j1j2・n_j1|＜R_j2sinθであれば、微小面要素j1を含む平面P1と微小面要素j2が交わっていると判定する。ここで、θは、２つの法線ベクトルn_j1とn_j2とがなす角度であり、θ=arcos(n_j1・n_j2)として表される。微小面要素j1を含む平面P1と微小面要素j2とが交わっていない場合（ステップＳ７３１否定）、反発力算出部３は、微小面要素j1と微小面要素j2とが交わらないと判定する（ステップＳ７３８）。

これに対して、微小面要素j1を含む平面P1と微小面要素j2とが交わっている場合（ステップＳ７３１肯定）、反発力算出部３は、微小面要素j2を含む平面P2と微小面要素j1が交わっているか否かを判定する（ステップＳ７３２）。微小面要素j2を含む平面P2と微小面要素j1とが交わっていない場合（ステップＳ７３２否定）、反発力算出部３は、微小面要素j1と微小面要素j2とが交わらないと判定する（ステップＳ７３８）。

これに対して、微小面要素j2を含む平面P2と微小面要素j1とが交わっている場合（ステップＳ７３２肯定）、反発力算出部３は、微小面要素j2の外縁と平面P1との交点を求める（ステップＳ７３３）。具体的には、反発力算出部３は、数式４で表される未知のベクトルrに対する連立方程式を解くことで、微小面要素j2の外縁と平面P1との交点３０１及び交点３０２を求めることができる。

数式４の連立方程式を解くと高々２つの交点座標が求まる。そこで、反発力算出部３は、求めた交点３０１及び交点３０２のうち少なくとも１つが微小面要素j1に含まれているか否かを判定する（ステップＳ７３４）。交点３０１及び交点３０２のうち少なくとも１つが微小面要素j1に含まれている場合（ステップＳ７３４肯定）、反発力算出部３は、微小面要素j1と微小面要素j2とが交わると判定する（ステップＳ７３７）。

これに対して、交点３０１及び交点３０２が１つも微小面要素j1に含まれていない場合（ステップＳ７３４否定）、反発力算出部３は、微小面要素j1の外縁と平面P2との交点を求める（ステップＳ７３５）。この場合、数式４のj1とj2を入れ替えた連立方程式を解くことで、反発力算出部３は、交点を求めることができる。

そして、反発力算出部３は、求めた交点のうち少なくとも１つが微小面要素j2に含まれているか否かを判定する（ステップＳ７３６）。交点のうち少なくとも１つが微小面要素j2に含まれている場合（ステップＳ７３６肯定）、反発力算出部３は、微小面要素j1と微小面要素j2とが交わると判定する（ステップＳ７３７）。

これに対して、交点が１つも微小面要素j2に含まれていない場合（ステップＳ７３６否定）、反発力算出部３は、微小面要素j1と微小面要素j2とが交わらないと判定する（ステップＳ７３８）。

以上のようにして、反発力算出部３は、図８のステップＳ７３における仮の最近接微小面要素と連続体粒子方向に突出するように交わる微小面要素があるか否かの判定を行うことができる。

［プログラム］
また、上記実施例で説明した各種の処理は、あらかじめ用意されたプログラムをコンピュータで実行することによって実現することができる。そこで、以下では、図１２を用いて、図２に示した連続体粒子の運動の解析処理を実行する機能を有する連続体運動解析プログラムを実行するコンピュータの一例を説明する。

図１２は、連続体運動解析プログラムを実行するコンピュータを示す図である。図１２に示すように、コンピュータ１０００は、ＲＡＭ(Random Access Memory)１０１０と、キャッシュ１０２０と、ＨＤＤ(Hard Disk Drive)１０３０と、ＲＯＭ(Read Only Memory)１０４０と、ＣＰＵ(Central Processing Unit)１０５０、バス１０６０とを有する。ＲＡＭ１０１０、キャッシュ１０２０、ＨＤＤ１０３０、ＲＯＭ１０４０、ＣＰＵ１０５０は、バス１０６０によって接続されている。

ＲＯＭ１０４０には、図２に示した連続体粒子の運動の解析処理を実行する機能を発揮する連続体運動解析プログラム１０４１が予め記憶されている。

そして、ＣＰＵ１０５０は、これらの連続体運動解析プログラム１０４１を読み出して実行する。これにより、図１２に示すように、連続体運動解析プログラム１０４１は、連続体運動解析プロセス１０５１になる。

なお、上記した連続体運動解析プログラム１０４１については、必ずしもＲＯＭ１０４０に記憶させなくてもよい。例えば、コンピュータ１０００に挿入されるフレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ(Compact Disk-Read Only Memory)、ＭＯ(Magneto-Optical)ディスク、ＤＶＤ(Digital Versatile Disk)ディスク、光磁気ディスク、ＩＣ(Integrated Circuit)カードなどの「可搬用の物理媒体」に連続体運動解析プログラム１０４１を記憶させてもよい。または、コンピュータ１０００の内外に備えられるハードディスクドライブ（ＨＤＤ）などの「固定用の物理媒体」に連続体運動解析プログラム１０４１を記憶させてもよい。または、公衆回線、インターネット、ＬＡＮ(Local Area Network)、ＷＡＮ(Wide Area Network)などを介してコンピュータ１０００に接続される「他のコンピュータ（またはサーバ）」に連続体運動解析プログラム１０４１を記憶させてもよい。そして、コンピュータ１０００は、上述したフレキシブルディスク等から各プログラムを読み出して実行するようにしてもよい。

以上に説明したように、本実施例では、固定境界を微小面要素で表し、各連続体粒子に対する微小面要素からの法線方向に向いた反発力を用いて固定境界から連続体粒子にかかる圧力が算出できる。微小面要素は適当なものを選択できるので、複雑な形状の固定境界を扱うことが容易となり、どのような形状の固定境界に対しても、本実施例に係る解析技術を適用することが可能となる。また、連続体粒子に対する微小面要素からの法線方向の反発力を用いて各連続体粒子に係る力を求めるので、単純な形状の固定境界であれば、連続体粒子に対して同じ方向に反発力が働くため、少数の微小面要素で反発力を正確に求められる。すなわち、少数の微小面境界で固定境界を表すことが可能となる。そのため、数値解析に必要なメモリ量と計算時間とをともに削減することが可能となり、同一の計算資源を用いた場合には、従来技術よりも大きなサイズの数値解析を行うことが可能となる。

また、本実施例では、反発力に対して各微小面要素と連続体粒子の影響範囲との交わる部分の面積に比例させた重み付けを行っている。これにより、微小面要素のサイズが一様でない形状の場合に、微小面要素同士の境目付近を連続体粒子が移動したとしても、反発力は連続的に変化するようになる。そのため、固定境界が連続体粒子に及ぼす影響をより現実的に表現することが可能となる。

また、境界面を微小円盤の集合として表現する場合、連続体粒子方向に突出する面を隙間無く覆うためには、一部の面要素が他の面要素と交わるという状況が必然的に生じる。この場合、連続体粒子方向に突出するように交わった面要素における境界の角の部分で反発力が不連続に変化するおそれがある。その際に、反発力が極めて大きな値を持つ場合があり、陽的な時間積分法を用いた数値解析では桁あふれ（オーバーフロー）が発生することにより、事実上計算をそれ以上進めることができなくなるという問題が従来技術にはあった。

この点、本実施例では、連続体粒子の受ける反発力のうち、無限大に発散する成分を最近接微小面要素からのみ受けるようにしている。これにより、連続体粒子方向に突出するように交わった面要素における境界の角付近で極端に大きな反発力を受ける機会を軽減し、陽的な解法を用いた場合にも数値計算を安定に進めることができる。

さらに、従来技術では、複数の境界をポリゴンで表したものを個別に用意して、それらを結合することによって連続体粒子に対する固定境界からの作用を表す１つの境界条件を設定する場合、複数の境界のつなぎ目の条件を適切に設定する必要がある。これに対して、本実施例では、境界面要素のつながりの計算は必要でないため、つなぎ目の条件を設定する必要が無くなり、複雑な境界条件を設定する手間を低減することができる。

以上のように、本願の連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置は、粒子法を用いた数値計算に有用である。特に、河川や海の水の流体解析に対して本願の連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置を用いることで防災の計画立案に役立てることができる。また、本願の連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置を用いて、鋳造過程を解析することで製品設計に寄与することができる。また、弾性体に対して、本願の連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置を用いることで、製品設計の際に封止ゲルの形状などを適切に決定することができる。

１ ユーザインタフェース
２ 記憶部
３ 反発力算出部
４ 圧力算出部
５ 解析部
６ 出力部

Claims (7)

1. 連続体を粒子として表した場合の各連続体粒子を表す影響範囲の情報を含むデータ、及び固定境界を任意の形状を有する微小な領域の集合として表した場合の法線方向の情報を含む、各微小面要素を表すデータの入力を受け付けるデータ受付手順と、
   前記連続体粒子の影響範囲と交わる前記微小面要素のそれぞれから前記連続体粒子に対して働く、前記微小面要素の法線方向の反発力を算出する反発力算出手順と、
   各連続体粒子に対してそれぞれの微小面要素から働く反発力を足し合わせることで、各前記連続体粒子に対して働く前記固定境界からの力を求める圧力算出手順と、
   各前記連続体粒子に対して働く前記固定境界からの力を基に、前記連続体の運動の解析を行う解析手順と
   をコンピュータに実行させる連続体運動解析プログラム。
2. 前記反発力算出手順は、前記連続体粒子と該連続体粒子の影響範囲に含まれる前記微小面要素との間の距離が少なくとも所定距離以上であれば、該距離に反比例する該微小面要素からの力を求め、さらに、求めた力に対して該連続体粒子の影響範囲と該影響範囲に含まれる前記微小面要素との交わった部分の面積に応じた重み付けを行うことで、該微小面要素から該連続体粒子に対して働く反発力を算出する請求項１に記載の連続体運動解析プログラム。
3. 前記反発力算出手順は、前記反発力の対象となる連続体粒子の影響範囲と交わる前記微小面要素のうちの少なくとも１つの微小面要素において、該連続体粒子と該微小面要素との距離が０になる場合には該微小面要素からの反発力が発散するように該微小面要素からの反発力を求めることを特徴とする請求項１又は請求項２に記載の連続体運動解析プログラム。
4. 前記反発力算出手順は、
   前記反発力の対象となる連続体粒子の影響範囲と交わる前記微小面要素を微小面要素集合とする集合決定手順と、
   微小面要素集合のうち、反発力の対象となる連続体粒子がその法線方向に位置し、且つ、連続体粒子との距離が最も近い微小面要素を仮の最近接微小面要素とする仮最近接決定手順と、
   該仮の最近接微小面要素が微小面要素集合に含まれる他の微小面要素と前記連続体粒子方向に突出するように交わる場合、該仮の最近接微小面要素とされた微小面要素を微小面要素集合から除外する除外手順と、
   前記仮最近接決定手順及び前記除外手順を繰り返し、微小面要素集合に含まれる他の微小面要素と前記連続体粒子方向に突出するように交わらない仮の最近接微小面要素を求め、求めた仮の最近接微小面要素を最近接微小面要素として特定する最近接特定手順と、
   該連続体粒子と該最近接微小面要素との距離が０になる場合には該最近接微小面要素からの反発力が発散するように該最近接微小面要素からの反発力を求める最近接反発力決定手順と
   を含むことを特徴とする請求項３に記載の連続体運動解析プログラム。
5. 前記解析手順は、特定の時刻における各前記連続体粒子に対して働く前記固定境界からの力を基に、各前記連続体粒子の加速度を求め、求めた加速度を時間で積分することで前記特定の時刻における前記連続体の動きを求め、さらに、各時刻における前記連続体の動きを時系列に沿ってまとめることで前記連続体の運動の解析を行うことを特徴とする請求項１に記載の連続体運動解析プログラム。
6. 連続体を粒子として表した場合の各連続体粒子を表す影響範囲の情報を含むデータ、及び固定境界を任意の形状を有する微小な領域の集合として表した場合の各微小面要素を表す法線方向の情報を含むデータの入力を受け付けるステップと、
   各前記連続体粒子の前記影響範囲と交わる前記微小面要素のそれぞれから各前記連続体粒子に対して働く前記微小面要素の前記法線方向の反発力を算出するステップと、
   各前記反発力を足し合わせることで各前記連続体粒子に対して働く前記固定境界からの力を求めるステップと、
   各前記連続体粒子に対して働く前記固定境界からの力を基に、前記連続体の運動の解析を行うステップと
   を有することを特徴とする連続体運動解析方法。
7. 連続体を粒子として表した場合の各連続体粒子を表す影響範囲の情報を含むデータ、及び固定境界を任意の形状を有する微小な領域の集合として表した場合の各微小面要素を表す法線方向の情報を含むデータの入力を受けて各データを記憶する記憶部と、
   各前記連続体粒子の前記影響範囲と交わる前記微小面要素のそれぞれから各前記連続体粒子に対して働く前記微小面要素の前記法線方向の反発力を算出する反発力算出部と、
   各前記反発力を足し合わせることで各前記連続体粒子に対して働く前記固定境界からの圧力を求める圧力算出部と、
   各前記連続体粒子に対して働く前記固定境界からの力を基に、前記連続体の動きの解析を行う解析部と
   を備えたことを特徴とする連続体運動解析装置。
   

Images