JP4253822B2

JP4253822B2 - 粒子法の解析プログラム

Info

Publication number: JP4253822B2
Application number: JP2007211559A
Authority: JP
Inventors: 昭夫堀; 拡美向上
Original assignee: 株式会社テクニカルスルー
Priority date: 2007-07-18
Filing date: 2007-07-18
Publication date: 2009-04-15
Anticipated expiration: 2027-07-18
Also published as: US7672821B2; US20090024376A1; JP2009026279A

Description

本発明は，物体を粒子で置換して，粒子が互いに接触していれば，粒子間の力や熱の伝達や物質移動を考慮して，粒子の変位や変形を計算する広義の粒子法のプログラムに関する。

有限要素法は，連続体を節点で囲まれた区分的な要素に分割して，要素間の主として物理的な関係を，節点間の関係で記述する解析方法である。このため有限要素法で接触問題を解析するためには工夫が必要で，例えば，重なり量を修正する力を計算させる手順を用いる特許文献１や，接触面間に新たな要素を発生させる特許文献２や，２要素が衝突する場合の衝突面の方向を両要素表面の平均方向で与える特許文献３がある。

一方，粒子法は，離散系の解析に本質的に向いているものの，現状では粒子自体に形状の概念がないか又は形状を球にしている（例えば特許文献４や非特許文献１や非特許文献２参照）。

特許２６９２６６８号特許３２４２８１１号特開２００２−２２９９６９号公報特開平０７−３３４４８４号公報越塚誠一著，「粒子法」，丸善，２００５年２月発行伯野元彦著，「破壊のシミュレーション」，森北出版，１９９７年１０月発行

以上に述べた従来の解析方法は，有限要素法を用いたものでは，連続体を基本に据えて一部要素の接触を扱うために離散系の解法には元々不向きになっている。また，例えば直方体の接触を考えればわかるように，接触形態には面と面だけでなく，面と角，辺と辺，それら様態間での移行，と種々の様態があるが，従来の有限要素法による方法では，これらの複雑な様態に対応できていない，という課題があった。

一方，粒子法を用いた方法では，粒子形状に球のような簡単な形状しか扱えず，より複雑な形状を模擬するには特に３次元中では多大な粒子数を要する，という問題があった。また多大な粒子数をもって模擬しても，微視的な形状に球形が残ったままでいると，粒子の集合体同士の接触方向が微視的な粒子形状で決定されてしまう，という課題があった。

これら課題を解決するための請求項１の発明は，粒子法において，粒子形状を球形でなくて構わなくするものである。その具体的手段は，各粒子に形状及び大きさを表す数学的に滑らかな閉曲面を設定する事と，近接する２粒子の接触の有無の判定を該閉曲面の接触の有無に基づいて行う事と，をまず含んでいる。数学的に滑らかな閉曲面にすると，基本的に２粒子間の接触が点になる可能性が非常に高く，接触の様態をまず単純化して考える事ができる。なお厳密には２粒子間の接触点数が２以上の可能性も残ってはいるが，本願発明では該接触点数が通常の場合すなわち１である場合に基づいて考える。

次に，球形でない粒子の移動量を計算するためには，該２粒子間に働く力の作用位置を定める必要がある。このため請求項１の発明では，該接触が有と判定された２粒子間の計算上の接触点を，該両閉曲面がちょうど１点で接するように該両閉曲面を同じ倍率で等方的に縮小した時の点として一意に定めている。もし本方法を用いない場合，例えば，以下のような問題が起きる。いま，任意の径長さを持つ２つの楕円体が，３次元中で任意の方向に置かれ，任意の相対位置関係で一部が重なった場合を考える。楕円体の表面である楕円面の重なりは，数式的に大変複雑な閉曲線となるから，以後の解析手法の構築が非常に困難となる。一方，本方法を用いれば，該両閉曲面の縮小率さえ数値的に求めておけば，縮小によって重なりを点にできるため，以後の解析手法の構築が非常に簡単になる。さらに，本方法の別の利点として，閉曲面形状を楕円面のような特別な形状に限定する必要がないので，粒子形状をより自由にできると共に，接触する２粒子の形状が全く異なっても構わなくなる。

球形でない粒子の移動量を計算するためには，該２粒子間に働く力の方向も定める必要がある。このため請求項１の発明では，該縮小両閉曲面であれば点接触する事，つまり，該計算上接触点では該縮小両閉曲面の法線ベクトルの方向が一致する事，通常は１８０度逆向きであるが，を利用して，該計算上接触点を通る平面で，かつ該法線方向を垂線方向とする平面，を計算上の接触平面として定めている。

該計算上接触点や該計算上接触平面を用いると，計算上の接触面積を容易に設定する事ができる。例えば，粒子の大きさから粒子体積を定めておき，粒子体積の半分を粒子中心点から該計算上接触平面までの距離Ｌで除した値を計算上の接触面積の候補値とし，該２粒子のそれぞれについて求めた候補値のうち小さい値を採用する方法で計算する事ができる。

該計算上接触点や該計算上接触平面や該計算上接触面積を用いると，該２粒子間の剛性や伝達抵抗や輸送抵抗又はそれらの逆数としての粒子間特性値又は特性値行列を計算する事が容易になる。例えば，第３の粒子がない場合に，前記距離Ｌの区間の剛性を，該計算上接触面積と粒子の材料剛性との積を前記距離Ｌで除した値で与えて，該２粒子のそれぞれの側の区間剛性に基づいて該２粒子間の軸剛性を計算する方法がある。

該粒子間特性値又は特性値行列をもとに，複数粒子の自由度間の特性値行列を計算すれば，有限要素法と類似の手続きにより，連立１次方程式を解いて，該複数粒子の並進や回転の移動量又は移動量増分を計算する事が可能になる。以上述べたように，請求項１の発明は前記課題を解決する事ができる。

請求項２の発明は，請求項１のプログラムを有限要素法のプログラムと連成させて構成した場合である。請求項１の発明に加えて，粒子の並進や回転の自由度又は自由度の一部を有限要素法の節点自由度と等置して，複数粒子の自由度間の特性値行列の全部又は部分を有限要素法で計算される広義の剛性行列の全部又は部分に加算する事で，請求項２の発明は有限要素法と簡単に共存可能で連成解析を可能にする。請求項２の発明は，本願の粒子法が有限要素法と二者択一的なものではなく，解析対象の連続体部分には有限要素を使う事で，より適用範囲の広い計算やより精確な計算ができ，過去のプログラム資産も本発明に有効利用できる事を示している。

請求項３の発明は，請求項１又は請求項２の発明において，該両閉曲面の縮小率や該計算上接触点を数値的に求める手順をより具体的に限定して構成した場合である。該計算上接触点を求める前にまず求めるべき該縮小率は単一のスカラーに過ぎないが，滑らかな閉曲面が点接触するという条件自体が，そもそも多変数の高次方程式に何らかの判別式を要求するという事を意味するため，該縮小率の計算手順には工夫が必要となる。請求項３は，その１例を示しており，以下の手順を用いている。最初に，計算上接触点の仮定点を仮定する。次に，該仮定点を通るように該両閉曲面を異なる倍率で等方的に拡大縮小する。該仮定点における該拡大縮小両閉曲面の両法線ベクトル方向及び両局部曲率を計算する。該両局部曲率を用いれば，該仮定点の過大でない位置変化に対する該両法線ベクトル方向の相違の変化率が求まるので，該仮定点位置を修正しながら数回の収束計算を行えば，該両法線ベクトル方向の相違を簡単に許容値以下にする事ができる。また，該仮定点の位置変化に対する該両拡大縮小率の相違の変化率も簡単に計算できるから，該収束計算中に同時に該両拡大縮小率の相違も許容値以下に収束させる事ができる。該収束計算によれば，該仮定点位置を該計算上接触点の位置に収束でき，該縮小率や該法線ベクトル方向も同時に得る事ができる。請求項３の発明は，以上のように該計算上接触点を数値的に一意に定める事ができる。

請求項４の発明は，請求項１〜請求項３の発明のいずれかにおいて，接触の有無の判定方法に請求項３の収束計算手順を応用して構成した場合である。前述と同様に多変数の高次方程式から何らかの判別式を求めて解の有無を調べるのは簡単でない場合が多いので，請求項３の収束計算手順を使って数値的に求めた該拡大縮小率に基づいて接触の有無を判定する事は，接触の有無の判定が簡単でない場合には特に有効である。該２粒子間で引張力を伝達しない場合には，該拡大縮小率が１未満の場合にのみ接触を考えれば良く，引張力を伝達する場合には略１以下で接触を考えれば良い。請求項４の発明は，接触の有無をより正しく数値的に判定する手順を含んでいる。

請求項５の発明は，請求項１〜請求項４の発明のいずれかにおいて，該閉曲面の大きさを表現するものとして，立体３方向の径又は径変化量又は径変化率を採用し，それを粒子の自由度の一部として構成した場合である。粒子の径変化を自由度に組み入れる事で，粒子の大きな寸法変化及びその影響を解析全体の中に直接的に取り込む事ができる。一般に，有限要素法は，節点間の変位差から要素内の歪が決まるので，小さな歪と大きな歪を特に区別する必要がないが，粒子法で粒子に大きさを設定すると，粒子の接触による重なりと粒子の大きさや形状の変化とは異なる意味あいを持つようになる。請求項５の発明は，粒子の径を変化させる事で，流体の流動や固体の塑性変形のような大きな変形を模擬する事を可能にしている。このような粒子の径変化は，請求項１のように粒子形状を球形でなくて構わなくした後に初めて可能にできる。また逆に，粒子の自由度である径変化を単に拘束するだけで平面歪のような状態を模擬できるため，より便利にもなる。

請求項６の発明は，請求項５の発明において，複数粒子の持つ自由度ベクトルに特性値行列を数式的又は数値的に前から乗じたベクトル，すなわち有限要素法の広義の外力ベクトルと数式的な位置関係が略対応するベクトル，が前記径や径変化量や径変化率の自由度と対応する行の要素に，偏差応力に相当する物理量を持つように構成した場合である。連立１次方程式を解いて粒子の変位を求めるには，広義の外力ベクトルの側が既知である必要があるが，粒子の径変化自由度と同じ行にある外力ベクトル側の値をどのような物理量にするかは解析方法の留意点の１つになる。しかしながら，例えば流体の流動の一部を粒子の形状変化を使って模擬する場合，該物理量を偏差応力にすると模擬しやすい。また固体の塑性変形も，一般に偏差応力によるとされている。以上の理由により，請求項６の発明は，該物理量を応力の１種である偏差応力として，径変化による各種現象の模擬をより容易にしている。

請求項７の発明は，請求項１〜請求項６の発明のいずれかにおいて，粒子の歪を計算する手順をやや限定して構成した場合である。一般に有限要素法では，要素の独立性が強く，要素節点の変位だけから歪を計算できるが，粒子法では，ある１つの粒子に接触する粒子の数自体が変化するので，該２粒子の接触だけから粒子の歪を計算する事ができない。このため請求項７の発明では，該縮小率に基づいて該接触による該粒子の歪への寄与量を求めると共に，該粒子に接触する全粒子からの該寄与量を該粒子についてスカラー合計又はベクトル合計して該粒子の歪を計算する手順を加えている。これにより，接触による粒子の歪をより適切に計算できる。

請求項８の発明は，請求項７の発明において，粒子間の力を計算する手順をやや限定して構成した場合である。粒子間の力を求める事は解析の基本的事項であるが，該計算上接触点及び該法線ベクトル方向及び該計算上接触面積及び前記歪及び材料の機械的性質を用いれば，粒子間の接触力の全部又は部分を計算する事が簡単にできる。

請求項９の発明は，以上では立体空間を主眼に説明してきた請求項１〜請求項８の発明が，基本的には平面解析にもそのまま適用できる事から，面を線に読み替えて次元数を１つ減らしたものである。なお読み替えは，閉曲面を閉曲線に，接触平面を接触線にする程度で良い。

本発明によれば，有限要素法が不向きである接触や離間を繰り返すような離散系の問題を数値解析できる。また，粒子形状が滑らかな閉曲面であれば良いため，接触の様態をより正しく模擬できると共に，少ない粒子数で複雑な形状の問題を解析できる。さらに，請求項２の発明によれば連続体部分に有限要素を用いて本粒子法と共存させた解析が，請求項５の発明によれば大変形を伴う解析が，それぞれできる。

以下，本発明の実施の最良の形態を図１〜図３に基づいて説明する。図１は最良の形態のフローチャートで，プログラム全体の流れを示す。以下，流れに沿って説明する。最初に入力データの取得を行う（Ｓ１）。該取得は，入力装置から行っても良く，既存のファイルから行っても良い。該入力データは，各粒子について，例えば，中心の座標，３方向の径，径の方向，形状を示す閉曲面の関数形，初速度，材料物性，温度，電荷，等で良い。但し該閉曲面の形状は，数学的に滑らかな形状とする。大きな変形も計算できるように，粒子の径変化率を解析全体の自由度の一部にする（Ｓ１）。広義の外力ベクトルにおける，該径変化率の自由度に対応する要素の物理量は偏差応力とする。

次に，近接する２粒子に関する接触の有無を，該閉曲面の接触の有無に基づいて判定する（Ｓ２）。該接触が有と判定されれば，以下の計算（Ｓ３〜Ｓ８）を行う。但し該判定が簡易な計算式から決まらない場合には，下記の計算上の接触点の計算（Ｓ３）の中で該判定を併せて行う。

計算上の接触点の計算（Ｓ３）は以下の収束計算により行う。図２は接触の様態の概念図，図３は計算上の接触点の計算（Ｓ３）のフローチャートである。最初に計算上の接触点３の位置を仮定する（Ｓ３１）。該仮定点３を通るように粒子Ａの閉曲面１ａ及び粒子Ｂの閉曲面１ｂを異なる倍率で等方的に拡大又は縮小する（Ｓ３２）。該両粒子Ａ及びＢの該拡大縮小両閉曲面１Ａ及び１Ｂについて，該仮定点３における両法線ベクトル方向５Ａ及び５Ｂ，及び両局部曲率を求める（Ｓ３３）。該両法線ベクトル方向５Ａ及び５Ｂは，該収束計算の途中では必ずしも平行ではないから，両者の相違を修正するように該仮定点３の位置を修正する（Ｓ３４）。必要な修正量は，該局部曲率が該法線ベクトル方向を該仮定点３の移動量で偏微分する意味を持つ事から，容易に求める事ができる。同時に，両閉曲面の該拡大縮小率の相違も修正するように，該仮定点３の位置を主に法線方向に修正する（Ｓ３４）。該両法線ベクトル方向５Ａ及び５Ｂの相違，及び該両拡大縮小率の相違が許容値以下になるまで，ステップＳ３２〜Ｓ３４を繰り返し行って収束計算する（Ｓ３５）。収束後には，該仮定点３が計算上の接触点３となる。なおステップＳ２の接触の有無の判定がスッテプＳ３に持ち越されている場合には，収束後の該拡大縮小率が１未満又は略１以下であれば接触有りと判定する。

次に，計算上の接触平面の計算（Ｓ４）を行う。計算上の接触平面４は，該計算上接触点３を通る平面で，収束後の該両法線ベクトル方向５Ａ及び５Ｂが垂線方向となるよう定める。

続いて，以下の方法で計算上の接触面積の計算（Ｓ５）を行う。まず，両粒子Ａ及びＢの中心点２Ａ及び２Ｂから該接触平面４までの距離５Ａ及び５Ｂをそれぞれ求める。次に，両粒子Ａ及びＢの片方から見た接触面積ＳＡ及びＳＢを，両粒子Ａ及びＢの持つ体積の半分ＶＡ及びＶＢをそれぞれ該距離５Ａ及び５Ｂで除した値で与える。計算上の接触面積は，ＳＡとＳＢのうちの小さい方の面積とする。小さい方にする理由は，大きい方の残された面積には，他の粒子が接触する可能性が十分あるからである。

次に，体積歪への寄与量の計算（Ｓ６）を行う。基本的には，ステップＳ３で得られた収束後の該拡大縮小率に基づき計算するが，両粒子Ａ及びＢの材料剛性の相違や，前記面積ＳＡとＳＢの相違の影響を考慮して計算しても良い。次に，体積歪への影響係数ベクトルの計算（Ｓ７）を行う。該影響係数ベクトルは，体積歪への寄与量を粒子の並進や回転や径変化率の各自由度で偏微分したものとし，計算過程はやや複雑になるものの，該計算上接触点３や，該法線ベクトル方向及び該距離の５Ａ及び５Ｂや，該面積ＳＡ及びＳＢを求める計算式において，該自由度に関する偏微分係数を順に数値計算してゆくと求める事ができる。

以上のステップＳ３〜Ｓ７の主要な計算結果を内部メモリー等の記憶装置に保存する（Ｓ８）。接触が有と判定された全ての２粒子の組合せについてステップＳ３〜Ｓ８を繰り返す。

ステップＳ６で求めた体積歪への寄与量を記憶装置から読み出して，該粒子が関係する全接触で合計して，該粒子の歪を求める。合計の方法は，流体ではスカラー合計，固体ではベクトル合計が良いが，解析目的に応じて選択可能とする（Ｓ９）。

ステップＳ３〜Ｓ５で求めた該計算上接触点及び，該法線ベクトル方向５Ａ及び５Ｂ，及び該計算上接触面積，及びステップＳ９で求めた該粒子の歪，及び粒子材料の剛性，を用いて粒子間の接触力を計算する（Ｓ１０）。なお粒子間の接触力には摩擦力などを別に加えても良く，また力だけでなく，粒子間の熱や物質や電荷の移動を計算しても構わない。該計算は接触する２粒子の組合せの数だけ行う。

有限要素があれば，要素の応力等を計算する。必要があれば，粒子や要素の慣性力や電磁気力も計算する（Ｓ１１）。

次に，ステップＳ７で求めた体積歪への影響係数ベクトルを記憶装置から読み出して，ステップＳ１０で用いた粒子間の接触力を求める式を使って，粒子間の剛性等を計算して広義の剛性マトリクスに足し込む。また，該接触力を求める式のうちの体積歪以外の部分を粒子の並進や回転や径変化率の各自由度で偏微分したものを使って，該偏微分に関わる粒子間の剛性等も計算して同様に広義の剛性マトリクスに足し込む（Ｓ１２）。なお広義の剛性マトリクスには，摩擦力や粒子間の熱や物質や電荷の移動に関わる伝達抵抗や輸送抵抗を別に加えても良い。

有限要素があれば，要素の剛性等を計算して広義の剛性マトリクスに足し込む。必要に応じて，粒子や要素の遠隔力に関する剛性等を計算して足し込む（Ｓ１３）。

広義の外力ベクトルを計算し（Ｓ１４），連立１次方程式を解いて，粒子の並進や回転や径変化率の各自由度や要素節点自由度の増分解を求める（Ｓ１５）。変位や応力等の計算結果のデータを記録装置のファイル等に出力する（Ｓ１６）。なお出力は別の時点で行っても良い。増分解析するため，以上のステップＳ２〜Ｓ１６を繰り返して計算実行する。該出力ファイルは，別に用意されたプログラムで必要な数値や図を画面に表示したり，紙に印刷したりする時に利用される。

本発明は，有限要素法を利用する産業分野において幅広く利用可能で，機械分野や土木建築分野や電機分野での製品設計や構築物設計のためのプログラムとして広範に利用可能である。また本発明を利用したプログラムそのものを市販可能である。

最良の形態のフローチャート接触の様態の概念図計算上の接触点の計算のフローチャート

符号の説明

１Ａ，１Ｂ粒子Ａ，Ｂの拡大縮小閉曲面
１ａ，１ｂ粒子Ａ，Ｂの閉曲面
２Ａ，２Ｂ粒子Ａ，Ｂの中心点
３計算上の接触点（収束するまでは仮定点）
４計算上の接触平面
５Ａ，５Ｂ粒子Ａ，Ｂの中心から計算上の接触平面までの距離（計算上の接触点を通る法線ベクトルでもある）

Claims

物理的や化学的な現象を粒子又は有限要素の集まりで表現し，球形でない粒子を１つ以上含む複数の粒子又は粒子及び要素節点の並進や回転の移動量を計算するプログラムであって，（ａ）近接する２粒子の接触の有無を，各粒子の形状及び大きさを表す数学的に滑らかな閉曲面の接触の有無に基づいて判定する手順と，該接触が有と判定された２粒子間に関する以下の（ｂ）〜（ｅ）の手順，（ｂ）該両閉曲面がちょうど１点で接するように該両閉曲面を同じ倍率で等方的に縮小した時の接触点を計算上の接触点として一意に定める手順，（ｃ）該計算上接触点を通る平面で，かつ該計算上接触点における該縮小閉曲面の法線ベクトルの方向を垂線方向とする平面を計算上の接触平面として求める手順，（ｄ）該計算上接触点や該計算上接触平面を用いて計算上の接触面積を定める手順，（ｅ）該計算上接触点や該計算上接触平面や該計算上接触面積を用いて該接触粒子間の剛性や伝達抵抗や輸送抵抗又はそれらの逆数としての粒子間特性値又は特性値行列を計算する手順と，（ｆ）該粒子間特性値又は特性値行列に基づいて複数粒子の自由度間の特性値行列を計算して，連立１次方程式を解いて該複数粒子又は該複数粒子及び要素節点の並進や回転の移動量又は移動量増分を計算する手順と，をコンピュータに実行させるためのプログラム。
前記手順ｆにおいて，粒子の並進や回転の自由度又は自由度の一部を有限要素法の節点自由度と等置する事により，複数粒子の自由度間の特性値行列の全部又は部分を，通常の有限要素法で計算される広義の剛性行列の全部又は部分に加算する手順，が挿入された請求項１のプログラム。
前記手順ｂの詳細が，（ｇ）該計算上接触点の仮定点を仮定する手順と，（ｈ）該仮定点を通るように該両閉曲面を異なる倍率で等方的に拡大縮小して，該仮定点における該拡大縮小両閉曲面の両法線ベクトル方向及び両局部曲率を求める手順と，（ｉ）該両法線ベクトル方向の相違及び該両拡大縮小率の相違を修正するように，該両局部曲率を用いて該仮定点の位置を修正する手順と，（ｊ）該両法線ベクトル方向の相違及び該両拡大縮小率の相違が許容値以下になるまで手順ｈ及び手順ｉを繰り返す収束計算の手順と，を含む請求項１又は請求項２のプログラム。
前記手順ａの詳細が，接触の有無の判定方法の１つとして，前記手順ｇ〜ｊを用いて求めた収束計算後の拡大縮小率が１未満又は略１以下の場合にのみ接触有りと判定する手順，を含む請求項１〜請求項３のいずれかのプログラム。
前記閉曲面の大きさを表現する立体３方向の径又は径変化量又は径変化率の全部又は一部が粒子の自由度の一部である，請求項１〜請求項４のいずれかのプログラム。
複数粒子の持つ自由度ベクトルに特性値行列を数式的又は数値的に前から乗じたベクトル，すなわち有限要素法の広義の外力ベクトルと数式的な位置関係が略対応するベクトル，が前記径や径変化量や径変化率の自由度と対応する行の要素に，偏差応力に相当する物理量を持つ，請求項５のプログラム。
前記手順ｂで適用した縮小率に基づいて該接触による該粒子の歪への寄与量を定める手順と，該粒子に接触する全粒子からの該寄与量を該粒子についてスカラー合計又はベクトル合計して該粒子の歪を求める手順と，をも実行させるための請求項１〜請求項６のいずれかのプログラム。
粒子間の接触力の全部又は部分を，粒子の前記歪及び粒子間の前記接触面積及び前記法線ベクトル方向及び粒子材料の機械的性質に基づいて計算する手順，をも実行させるための請求項７のプログラム。
以上の記載に対して，平面解析に適用可能なよう面を線に読み替えて次元数を１つ減らした，請求項１〜請求項８のいずれかのプログラム。