JP2007293540A - 連続弾性体変形シミュレーション方法、プログラム、および記録媒体 - Google Patents

Abstract

【課題】離散要素近似力学系モデルにおいて、誤差の少ない連続弾性体変形シミュレーション方法を提供する。
【解決手段】連続弾性体を隣接している８つの離散質量点をそれぞれ頂点とする６面体を単位キューブに分割し（８）、各単位キューブに対し、連続体としての直交する３軸の各軸方向の荷重に対する連続体面変位を計算し（１８）、各単位キューブに対し、各軸方向の荷重ベクトルに対する各頂点の変位を剛性マトリクスを用いて計算し（１４）、各単位キューブごとに、同一面内の上記頂点変位を平均して近似面変位を求め（２２）、各単位キューブごとに対応軸面の連続体面変位と近似面変位との差の絶対値が予め決められた誤差以下になるように剛性マトリクスの要素を調整する（２６）。
【選択図】図１

Description

この発明は例えば、人間や動物の皮膚軟組織つまり、舌や顔、あるいは木材などの建造物などの有体物、つまり弾性を有する連続体に荷重をかけた際、荷重を掛けられた部分がどのくらい変位するかをシミュレーションする連続弾性体変形シミュレーション方法、プログラムおよび記録媒体である。

人間は舌、顎、唇などの複数の調音器官や顔の筋をたくみに操ることで、調音運動を行っている。このような調音運動を調べるため、電気的、磁気的、光学的検出装置を用いた調音器官の代表点の位置計測や調音器官の運動生成に関与する筋の筋電図信号の計測やＭＲＩ（磁気共鳴画像）による形状計測などが行われてきた。近年ＭＲＩ動画撮影の技術も開発され、調音器官全体をカバーする時間・空間的な計測も可能になりつつある。
しかし、骨格系と異なり、弾性体である舌や唇の三次元運動は複雑に配置された筋によって組み立てられており、限られた数の代表点や筋電図信号で理解することは難しく、またＭＲＩ動画も細やかな筋の動きを観察するためには空間・時間分解能がまだ十分ではない。

このような観測的アプローチに対し、調音器官の形状や筋の配置および筋ダイナミックスを模擬する生理学的モデルを用いて、調音器官で生成される巧みな運動がどのように実現されるかを調べる構成的アプローチも検討されている。このアプローチでは、計測困難な複雑に入り込んだ筋の活動と運動との関係を陽に知ることが出来るという効果がある。

そしてＦＥＭ（有限要素法）により、正中矢状断面の舌の２次元モデルが構築されている。そして同様の手法を用いて、舌の３次元モデルも構築されている。また、ＦＥＭを用いて軟組織をモデル化するために、精密な数学的方法が提案され、それに基づいて舌の３次元モデルの構築が試みられている。しかし、これらＦＥＭを用いたモデルは、舌の変形を記述することができるものの、計算量が多く、形状が複雑な場合は、現実的でないという問題が指摘されていた。これらの方法に対し、複数の離散質量点に網状の弾性体である離散要素を張って、連続弾性体を構成する離散要素近似力学系モデルとして処理することが提案された。離散要素近似力学系モデルとは、非特許文献１〜４に記載のように、静力学／動力学的解析において、連続体の対象物を複数の離散質量点およびそれらを連結しあう離散要素で近似する手法であり、この手法により、モデリングの容易さや計算時間短縮などの効果が得られる。
Ｙ．Ｌｅｅ、Ｄ．Ｔｅｒｚｏｐｏｕｌｏｓ ａｎｄ Ｋ．Ｗａｔｅｒｓ，"Ｒｅａｌｉｓｔｉｃ ｍｏｄｅｌｉｎｇ ｆｏｒ ｆａｃｉａｌ ａｎｉｍａｔｉｏｎ，"Ｐｒｏｃ．ＳＩＧＧＲＡＰＨ９５，ＡＣＭ ＳＩＧＧＲＡＰＨ，ｐｐ．５５−６２，ＬｏｓＡｎｇｅｌｅｓ，ＣＡ，Ｕ．Ｓ．Ａ．，Ａｕｇ．１９９５． Ｒ．Ｇｒｚｅｓｚｃｚｕｋ ａｎｄ Ｄ．Ｔｅｒｚｏｐｏｕｌｏｓ，"Ａｕｔｏｍａｔｅｄ ｌｅａｒｎｉｍｇ ｏｆ ｍｕｓｃｌｅ−ａｃｔｕａｔｅｄ ｌｏｃｏｍｏｔｉｏｎ ｔｈｒｏｕｇｈ ｃｏｎｔｏｒｏｌ ａｂｓｔｒａｃｔｉｏｎ，"Ｐｒｏｃ．ＳＩＧＧＲＡＰＨ９５，ＡＣＭ ＳＩＧＧＲＡＰＨ，ｐｐ．６３−７０，ＬｏｓＡｎｇｅｌｅｓ，ＣＡ，Ｕ．Ｓ．Ａ．，Ａｕｇ．１９９５． Ｊ．Ｄａｎｇ ａｎｄ Ｋ．Ｈｏｎｄａ，"Ｃｏｎｓｔｒｕｃｔｉｏｎ ａｎｄ ｃｏｎｔｒｏｌ ｏｆ ａ ｐｈｙｓｉｏｌｏｇｉｃａｌ ａｒｔｉｃｕｌａｔｏｒｙ ｍｏｄｅｌ，"Ｊ．ｏｆ Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．，１１５（２），ｐｐ．８５３−８７０，Ｆｅｂ．２００４． 野添潤一，五味裕章，党建武，本多清志，"リアルな発話運動を実現する生理学的唇力学モデルの構築，"信学論（Ｄ−II），ｖｏｌ．Ｊ８８−Ｄ−II，ｎｏ．９，ｐｐ．１９４４−１９５３，Ｓｅｐ．２００５．

従来の上記離散要素近似力学系モデルでは、複数の質量点に網状の弾性体である離散要素を張って構成した近似力学系モデルにおいて、モデルを構成する離散要素の剛性として、非特許文献１〜３では一定値とし、非特許文献４では離散要素の長さを変数とする弾性係数を設定していた。しかし、連続弾性体をこのように近似した離散力学系は、離散要素の密度（単位空間体積に存在する離散要素の数）が一様でないと、静的荷重による質量点の変位は、実際の連続弾性体の変位と一致せず、誤差が生じてしまう。

この発明によれば、連続弾性体を複数の離散質量点に網状の離散要素を張った離散要素近似力学系モデルの連続弾性体変形シミュレーション方法において、隣接している８つの上記離散質量点をそれぞれ頂点とする６面体を単位キューブとし、各単位キューブに対し、キューブ内で直交する３軸の各軸方向の荷重に対する連続体面変位を計算し、一方、各単位キューブに対し、上記離散要素近似力学系における各軸方向の荷重ベクトルに対する各頂点の変位を剛性マトリクスを用いて計算し、この場合、上記荷重ベクトルとして、上記連続体面変位計算における対応する単位キューブの同一軸方向の荷重が各質量点に分割されて印加される構成とされ、各単位キューブごとに、同一面内の上記頂点変位を平均して近似面変位を求め、各単位キューブごとに対応軸面の上記連続体面変位と上記近似面変位との差の絶対値が予め決められた誤差以下になるように上記剛性マトリクスの要素を調整する。

以上の構成によれば、連続体面変位と近似面変位との差の絶対値が予め決められた誤差以下になるように剛性マトリクスの要素を調整するため、連続弾性体を近似している離散要素近似力学系モデルにおいて、モデリングの容易さや、計算時間短縮などの効果を維持しつつ、連続弾性体との適合性が成り立ち、正確な測定結果が得られる。

実施例１
この発明のハードウェア構成例を図１に示す。この実施例において、使用者がシミュレーションを所望する対象物を連続弾性体と定義する。入力手段２は連続弾性体入力手段４と荷重入力手段６とで構成されている。
まず使用者が連続弾性体を連続弾性体入力手段４により、入力し、荷重入力手段６に連続弾性体に与えたい荷重ｆ_ｃを入力する。連続弾性体入力手段４への入力は、例えば、ステレオカメラなどの、光学計測装置などで、対象である連続弾性体を読み取る手法や、図面がある場合、図面から直接読み取る手法などで、得た情報を数値化したデータとしてコンピュータに入力する。連続弾性体入力手段４は対象の連続弾性体の大きさ、形状が離散的な点として入力される。この離散的な点は従来の離散的質量点に網状の離散要素を張った離散要素近似力学系モデルにおけるノード、つまり離散質量点と対応している。入力された連続弾性体を示す情報、つまり各質量点を表す情報は、単位キューブ分割手段８に入力され、対象の連続弾性体が単位キューブに分割される。ここで、単位キューブとは、図２に示すようにｘ、ｙ、ｚ軸上で考える６面体であり、本発明では、この単位キューブを最小単位と考え、連続弾性体は単位キューブが集合したものとして考える。連続弾性体の形状が複雑な部分はその複雑さに比例して小さい単位キューブとされ、形状が単純な部分は大きな単位キューブとされる。つまり、各単位キューブはそれ自体が６面体の連続弾性体として、各直交３軸の各方向の荷重を受けた時に、その荷重に基づくその対向の弾性体面変位をフックの法則により、計算できる程度の大きさ、形状とする。また、連続弾性体を複数の離散質量点に網の離散要素を張った従来の離散要素近似力学系モデルにおける隣接している８個の質量点をそれぞれ頂点とする６面体は、離散要素近似力学系における単位キューブとなる。

分割された単位キューブが任意の順番に単位キューブ近似手段１０に入力される。入力された単位キューブは、図３に示すように、離散要素近似力学系モデルに近似される。なお、簡略化のため、以下は離散要素近似力学系モデルを離散モデルと言及する。
離散モデルにおける単位キューブはｘ、ｙ、ｚ軸上に質量点１〜８と、１つの質量点から他の７つの質量点とを繋ぎ合わされている２８本の離散要素から構成されており、各々の離散要素の剛性は弾性係数ｋである。離散モデルの構成情報は剛性マトリクス生成手段１２に入力されて離散モデルの剛性マトリクスＫが生成される。剛性マトリクス（Ｓｔｉｆｆｎｅｓｓ Ｍａｔｒｉｘ）とは、いわば、複数の弾性要素で構成した力学系において、単体の弾性要素における弾性係数と同じ役割を示すものであり、離散要素力学系を構成する弾性要素ｘの弾性係数ｋ^ｘ及びその方向性によって作られる。また剛性マトリクスＫは、２つ以上の要素（部材）を連結して、構築した構造物の外部荷重に対する挙動をマトリクス変位法もしくは剛性マトリクス法を用いて、解析する場合、フックの法則における変位の係数項としてでてくるマトリクスのことである。

なお、剛性マトリクス生成手段１２においての剛性マトリクスの生成方法に関しては、説明を簡略化にするため、２次元空間の場合について説明する。
図４に示すように離散要素力学系（トラス構造）は４個のノード（質量点）１〜４及びある１つのノードから他の３つのノードへ繋がっている６個の弾性要素（離散要素）Ａ〜Ｆで構成している。図４の力学系に対応する剛性マトリクスにおける、例えば弾性要素Ａの寄与分は、式(１)のように剛性マトリクスの中で展開される。以下の式(１)のマトリクスに付けておいた行番号及び列番号は、図４におけるノード番号に対応するものである。

ここで、式(１)のマトリクスの要素として展開されるｋ_０ ^Ａは、以下の式（２）のように定義する２×２（２行２列）のサブマトリクスである。

ただし、式（２）のｋ^ｍ（ｍ＝Ａ、．．．、Ｆ）は弾性要素Ａ〜Ｆの弾性係数を表し、ｅ_Ａ、ｆ_Ａは弾性要素Ａの方向余弦（ｄｉｒｅｃｔｉｏｎ ｃｏｓｉｎｅ）で、以下の式（３）〜（５）で表すことができる。
ｅ_Ａ＝（ｘ_２−ｘ_１）／ｒ_Ａ・・・（３）
ｆ_Ａ＝（ｙ_２−ｙ_１）／ｒ_Ａ・・・（４）
ｒ_Ａ＝｛（ｘ_２−ｘ_１）＋（ｙ_２−ｙ_１）｝^１／２・・・（５）
ここで、ｒ_Ａを弾性要素Ａの長さ、要素Ａの両端を構成するノード１及びノード２の座標をそれぞれ（ｘ_１、ｙ_１）、（ｘ_２、ｙ_２）とする。

次に弾性要素Ｂの寄与分は、弾性要素Ａの場合と同様、その接続関係から式（６）のように展開される。

剛性マトリクスの構成における個々の弾性要素の寄与分には、線形重ね合わせの原理が適用されるため、弾性要素Ａ及びＢの寄与分を考慮した剛性マトリクスは以下の式（７）のようになる。

弾性要素Ｃ〜Ｆに対しても剛性マトリクスの構成における寄与分を上記のように展開していくと、図４の離散要素力学系に対応する剛性マトリクスは以下の式（８）のように導出される。

式（８）の剛性マトリクスは式（２）のように定義する２×２のサブマトリクスｋ^Ａ _０〜ｋ^Ｆ _０の組み合わせを要素として持つため、その次元数は８×８＝６４次元（８行８列）となる。上記のように剛性マトリクス生成手段１２での処理を弾性要素の3次元配置を持つ離散要素力学系に対して拡張し、適用することによって３次元配置の剛性マトリクスＫを生成することが出来る。なお、３次元配置の剛性マトリクスの生成方法の詳細は、「構造工学の基礎と応用−例題で学ぶ 宮本裕 枝報堂出版 ２００３年５月」や「建設構造力学II 山田孝一郎他 森北出版」などに記載されている。

生成された剛性マトリクスＫは頂点変位計算手段１４に入力される。
頂点変位計算手段１４では、図３に示す離散モデルにおけるｘ軸、ｙ軸、ｚ軸方向それぞれの変位ベクトルδ^→ _ｘ、δ^→ _ｙ、δ^→ _ｚを求める。ここで、Ａ^→はベクトルＡを示すものとし、具体的計算方法を以下に示す。
荷重入力手段６よりの荷重ｆ_ｃが入力されると、成分分解手段２０でｘ軸、ｙ軸、ｚ軸方向の成分、ｆ_ｃｘ、ｆ_ｃｙ、ｆ_ｃｚに分解される。従って、ｘ軸方向成分の荷重ベクトルをｆ_ｘｘとする。ここで、添え字「ｘｘ」の、１つめのｘはｘ軸と垂直な断面（以下、ｘ軸断面という）を示し、２つめのｘはｘ軸方向を示すものとする。このｆ_ｘｘ ^→の各構成要素は次式で表せる。

ｆ_ｘｘ ^→＝［質量点１のｘ軸方向の荷重、質量点１のｙ軸方向の荷重、質量点１のｚ軸方向の荷重、質量点２のｘ軸方向の荷重、質量点２のｙ軸方向の荷重、質量点２のｚ軸方向の荷重、．．．、質量点８のｘ軸方向の荷重、質量点８のｙ軸方向の荷重、質量点８のｚ軸方向の荷重］
離散要素近似力学系では、図３中の質量点１〜４の変位を考える場合、質量点５〜８を固定されているものとし、各質量点１〜４は、荷重ｆ_ｃのｘ軸方向の成分ｆ_ｃｘを４等分したｘ軸方向の集中荷重が印加されるとする。つまり、ｆ_ｃｘは離散要素近似力学系において、以下の式（９）で表す１×１２のベクトルｆ_ｘｘ ^→と等価的に表すことが出来る。
ｆ_ｘｘ ^→＝［ｆ_ｃｘ／４、０、０、ｆ_ｃｘ／４、０、０、ｆ_ｃｘ／４、０、０、ｆ_ｃｘ／４、０、０］^Ｔ・・・（９）

また、質量点１〜４のｘ軸方向にかかる荷重をｆ_ｃｘ／４としているが、これは単位キューブを十分小さいものであると近似しているため、単位キューブが直方体や立方体に限らず、あらゆる６面体で適用される。
また離散要素近似力学系における荷重による弾性変形は、マトリクス変位法によって、求めることが知られている。マトリクス変位法の適用で、各質量点の変位を表すベクトルδ^→は剛性マトリクスＫの逆行列に各質量点に印加される変位力ベクトルｆとの積、つまり、以下の式（１０）のような線形連立方程式で求めることができる。
δ^→＝Ｋ^−１・ｆ^→・・・（１０）

ここで剛性マトリクスＫの簡略構成例を図５に示す。剛性マトリクスＫの構成要素ａ_ｍｎはｍ列ｎ行に位置することを示し、ａ_ｍｎは式（２）に対応する３×３のサブマトリクスである。よってこの剛性マトリクスＫは２４×２４の剛性マトリクスとなる。図３に示す例えば、ｘ軸方向の変位力（集中荷重）に対し、変位する面をｘ軸断面とすると、そのｘ軸断面の４つの各質量点は１、２、３、４、である。ある軸断面が質量点ｒ、ｓ、ｐ、ｑとして、構成される場合、その軸断面の変位を考える場合は、図５に示した剛性マトリクスにおいて、ｒ、ｓ、ｐ、ｑ行およびｒ、ｓ、ｐ、ｑ列以外の構成要素、を削除して、縮小した剛性マトリクス（Ｒｅｄｎｃｅｄ Ｓｔｉｆｆｎｅｓｓ Ｍａｔｒｉｘ）を考えればよい。従って、縮小された剛性マトリクスは１２×１２のマトリクスとなる。例えば、質量点１、２、３、４で構成された断面（ｘ軸断面）の変位（ｘ軸方向の変位）を考える場合、剛性マトリクスは図５Ｂに示すように構成要素ａ_ｍｎ（ｍ＝１、２、３、４、ｎ＝１、２、３、４）で構成される１２×１２の剛性マトリクスとなる。また、質量点２、３、５、８で構成された断面（ｙ軸断面）の変位（ｙ軸方向変位）を考える場合は図５Ｃに示すように構成要素ａ_ｍｎ（ｍ＝２、３、５、８、ｎ＝２、３、５、８）で構成される剛性マトリクスとなる。なお、上述したとおり、ａ_ｍｎは３×３のサブマトリクスである。

従って、例えば、図３中の質量点１、２、３、４で構成される断面をＢ_１２３４と表し、ｘ軸方向の断面Ｂ_１２３４の各質量点の変位を表す変位ベクトルδ_ｘ ^→ _１２３４を以下のように求めることが出来る。式（９）で示したｆ_ｘｘ ^→をｆ^→とし、このｆ^→と剛性マトリクス生成手段１２で生成された剛性マトリクスＫとを、式（１０）に代入すると、以下の式（１１）に示すようにｘ軸方向についての１２×１で表す質量点変位ベクトルδ_ｘ ^→ _１２３４を求めることが出来る。なお、このｘ軸方向の各質量点変位を求める例では、剛性マトリクスＫは図５Ｂに示したように縮小されたものである。
δ_ｘ ^→ _１２３４＝［δ_ｘ（１）、０、０、δ_ｘ（２）、０、０、δ_ｘ（３）、０、０、δ_ｘ（４）］^Ｔ・・・（１１）
なお、ベクトルの構成要素であるδ_ｉ（ｍ）は質量点ｍ（ｍ＝１、．．．、８）のｉ軸方向（ｉ＝ｘ、ｙ、ｚ）の変位を示し、Ａ^ＴはベクトルＡの転置行列を示す。

頂点変位計算手段１４で計算されたδ_ｘ ^→ _１２３４は頂点変位記憶手段１６に記憶される。以下、同様の処理での質量点２、３、５、８で構成される断面のｙ軸方向の質量点変位ベクトルδ_ｙ ^→ _２３５８、質量点３、４、７、８で構成される断面のｚ軸方向の質量点の変位ベクトル、δ_ｚ ^→ _３４７８が頂点変位計算手段１４で計算され、頂点変位記憶手段１６に記憶される。以下の式（１２）（１３）にδ_ｙ ^→ _２３５８、δ_ｚ ^→ _３４７８を示す。
δ_ｙ ^→ _２３５８＝［０、δ_ｙ（２）、０、０、δ_ｙ（３）、０、０、δ_ｙ（５）、０、０、δ_ｙ（８）０］^Ｔ・・・（１２）
δ^ｚ→ _３４７８＝［０、０、δ^ｚ _（３）、０、０、δ^ｚ _（４）、０、０、δ^ｚ _（７）、０、０、δ^ｚ _（８）］^Ｔ・・・（１３）

この発明の離散モデルにおいて、４個の質量点で構成される断面の変位δ_ｉ（以下、近似面変位といい、ｉ＝ｘ、ｙ、ｚとする）は当該４個の質量点の変位の平均値として求める。以下の式（１４）は図３の離散モデルの単位キューブに式（９）で示す荷重を加えた場合のｘ軸断面の近似面変位を求めるものである。
δ_ｘ＝１／４・Σ^４ _ｋ＝１δ_ｘ（ｋ）・・・（１４）
よって式（１１）、式（１４）からδ_ｘは式（１５）に示すように、
δ_ｘ＝１／４（δ_ｘ（１）＋δ_ｘ（２）＋δ_ｘ（３）＋δ_ｘ（４））・・・（１５）
となる。以下、同様の処理で同様にδ_ｙ、δ_ｚを求めることが出来る。そしてδ_ｘ、δ_ｙ、δ_ｚはそれぞれ近似面変位記憶手段２４に記憶される。

一方、成分分解手段２０によりｘ軸、ｙ軸、ｚ軸方向に分解された荷重は連続体面変位計算手段１８に入力される。また、単位キューブ分割手段８により分割された単位キューブが連続体面変位計算手段１８に入力される。連続体面変位計算手段１８では局所座標系においてある単位キューブのｉ軸断面にてｉ軸方向の荷重ｆ_ｃｉによるｉ軸方向の変位δ_ｃｉが以下の式（１６）により求められる。なお、式（１６）はフックの法則である材料力学の基本式として、広く知られ、ｉはｘ、ｙ、ｚを示し、つまり、ｘ、ｙ、ｚ軸方向それぞれの単位キューブの連続弾性体として変位が求められる。
δ_ｃｉ＝ｆ_ｃｉ・ｒ_ｃｉ／Ｅ・Ａ_ｃｉ・・・・（１６）

ここで、Ｅは連続弾性体の縦弾性係数（ヤング率）、ｒ_ｃｉは単位キューブのｉ軸方向の長さ、Ａ_ｃｉはｉ軸断面の面積を示す。ｘ、ｙ、ｚ軸方向それぞれの単位キューブの連続弾性体の変位（以下、連続面変位という）δ_ｃｘ、δ_ｃｙ、δ_ｃｚがそれぞれ求められると連続体面変位記憶手段１９に記憶される。

連続体面変位記憶手段１９よりの連続体面変位δ_ｃｘ、δ_ｃｙ、δ_ｃｚと近似面変位記憶手段２４よりの近似面変位δ_ｘ、δ_ｙ、δ_ｚがそれぞれ、調整手段２６に入力される。この発明での調整手段２６では、各軸方向における近似面変位の連続面変位に対する差の絶対値が、使用者が任意に設定した誤差εの範囲内になるように、つまり離散モデルは連続弾性体系と静的変形における幾何学的適合性をもつように剛性マトリクスＫの各要素、すなわち弾性係数を調整する。

具体的には以下の式（１７）（１８）を満たすよう、離散モデルを構成する上記の１２×１２の剛性マトリクスＫの各弾性要素の弾性係数を調整することによって、誤差εが小さい離散モデルが導出される。
│δ_ｃｘ−δ_ｘ│＋│δ_ｃｙ−δ_ｙ│＋│δ_ｃｚ−δ_ｚ│≦ε・・・（１７）
ただし、εは使用者が任意に設定した誤差の範囲である。
また式（１７）の代わりに、以下の式（１８）に示す、変位偏差の２乗を用いても、同様な効果をもたらす離散要素近似力学系が求められる。
ｃ_１（δ_ｃｘ−δ_ｘ）^２＋ｃ_２（δ_ｃｙ−δ_ｙ）^２＋ｃ_３（δ_ｃｚ−δ_ｚ）^２≦ε・・・（１８）

なお、ｃ_１、ｃ_２、ｃ_３は任意の定数であり、ｃ_１＝ｃ_２＝ｃ_３＝１／２にするのが好ましく、δ_ｃｘ、δ_ｃｙ、δ_ｃｚはそれぞれ固定値であることに留意されたい。式（１７）、あるいは式（１８）を満たすように、離散モデルにおける弾性要素の弾性係数の調整には、Ｎｅｗｔｏｎ−Ｒａｐｓｏｎ法、Ｂｉｓｅｃｔｉｏｎ法、最急降下法（ｓｔｅｅｐｅｓｔ ｄｅｓｃｅｎｔ ｍｅｔｈｏｄ）などあらゆる最適化手法の適用が可能である。また上記式（１５）（１７）（１８）は、＋ｘ、＋ｙ、＋ｚ軸断面に関して言及したが、−ｘ、−ｙ、−ｚ軸断面について適用しても問題はない。

式（１７）あるいは式（１８）で算出されるδ_ｘ、δ_ｙ、δ_ｚが出力手段２８から出力される。δ_ｘ、δ_ｙ、δ_ｚはそれぞれδ_ｃｘ、δ_ｃｙ、δ_ｃｚと近似するよう調整された値なので、連続弾性体の変位に対し誤差が小さい近似面変位をシミュレーションをすることができる。なお、式（１７）又は式（１８）を満足する調整が得られた時のδ_ｘ、δ_ｙ、δ_ｚがそれぞれ求める近似面変位である。

実施例２
この発明の実施例２では、また、実施例１で示したハードウェア構成に対し、破線で示すように共有要素平均手段３０が付加され、出力手段２８に修正値頂点変位計算手段３４が含まれる。その他の部分は実施例１と同様の処理を行う。
また、調整手段２６による調整終了後の離散モデルの弾性要素を連続体適合弾性要素として定義し、連続体適合弾性要素で構成した単位キューブを連続体適合単位キューブと定義し、連続体適合単位キューブの組み合わせで構成した全体系を連続体適合離散要素近似力学系として定義する。また調整手段２６により調整された剛性マトリクスを第１の剛性マトリクスＫ’と定義する。

第１の剛性マトリクスＫ’は、共有要素平均手段３０に入力される。共有要素平均手段３０では、連続体適合離散要素近似力学系にて、複数の単位キューブに共有される離散要素に対しては、各共有キューブから求めた連続体適合弾性係数の平均値をその要素の連続体適合弾性係数として設定する。具体的に、図６Ａ、Ｂを参照しながら説明すると、連続体適合単位キューブＡ、Ｂ、Ｃ、Ｄ（以下、適合単位キューブという）について適合単位キューブＡの右側面の前方下端の頂点（質量点）をａ１、後方下端の頂点をａ６、前方上端の頂点をａ１０、左側面の前方上端の頂点をａ５、前方下端の頂点をａ１１、後方上端の頂点をａ１２、適合単位キューブＢの右側面の前方上端の頂点をａ２、後方上端の頂点をａ７、適合単位キューブＣの左側面の前方下端の頂点をａ３、後方下端の頂点をａ８、適合単位キューブＤの左側面の前方上端の頂点をａ４、後方上端の頂点をａ９、とする。

適合単位キューブＡ、Ｂの各右側面と適合単位キューブＣ、Ｄの各左側面を互いに接合して、図６Ｂに示すような、適合単位キューブＡ、Ｂ、Ｃ、Ｄが組み合わされた形態を考える。
適合単位キューブＡのａ１とａ６とを結ぶ離散要素ａ１ａ６と、適合単位キューブＢのａ２とａ７とを結ぶ離散要素ａ２ａ７と適合単位キューブＣのａ３とａ８とを結ぶ離散要素ａ３ａ８と適合単位キューブＤのａ４とａ９とを結ぶ離散要素ａ４ａ９とは互いに共有関係にある。よって、離散要素ａ１ａ６の弾性係数ｋ_１６と離散要素ａ２ａ７の弾性係数ｋ_２７と離散要素ａ３ａ８の弾性係数ｋ_３８と離散要素ａ４ａ９の弾性係数ｋ_４９との平均値、つまり（ｋ_１６＋ｋ_２７＋ｋ_３８＋ｋ_４９）／４を取る。

また、例えば、適合単位キューブＡの質量点ａ５に接続される離散要素ａ５ａ１０、ａ５ａ１１、ａ５ａ１２、は他の適合単位キューブＢ、Ｃ、Ｄ中の離散要素のいずれとも共有されていないので、これらの弾性係数は変更しない。
これら、各共有離散要素に対し求めた修正弾性係数と、共有されずに変更されない弾性係数を用いて構成された剛性マトリクスを第２の剛性マトリクスＫ’’とすると、第２の剛性マトリクスＫ’’は連続体適合離散要素近似力学系の剛性マトリクスである。第２の剛性マトリクスＫ’’は修正値頂点変位手段３４に入力され、修正値頂点変位手段３４で、連続体適合離散要素近似力学系に対する頂点変位を計算する。

具体的には、汎用的に使用できる式（９）に、式（１０）で表すｆ_ｘｘ ^→とＫ’’を入力して、連続体適合単位キューブのｘ軸、ｙ軸、ｚ軸方向の変位であるδ_ｘ’、δ_ｙ’、δ_ｚ’を算出して出力する。共有離散要素の平均をとって第２の剛性マトリクスＫ’’を算出することで、より適合性のとれた、つまり、より誤差の少ない変位を算出することが出来る。

実験結果
単純な幾何学的形状を持つ連続弾性体に対して、調整手段２６による調整と共有要素平均手段３０による修正がされた連続体適合離散要素近似力学系を適用した結果を以下に説明する。図７Ａ、Ｂは立方体の連続弾性体を４個の同一寸法の単位キューブで構成した離散要素近似力学系に対する静的変形の結果であり、図７Ａは従来の技術の変形結果であり、図７Ｂはこの発明を適用した変形結果である。なお、両者とも連続弾性体の寸法は０．１×０．１×０．１（ｍ^３）で、その縦弾性係数（ヤング率）は１００（ｋＰａ）として設定した。荷重条件としては、+ｘ軸断面上における１２（ｋＰａ）の一様引張り応力分布を適用した。この発明の適用に当たっては誤差εを１０^−４（ｍ）とした。直線で示されているのは離散モデルでの単位キューブの変形を示し、２点鎖線で示されているのは、連続弾性体の変形を示している。

図７Ａにおいて、一様応力分布のため、連続弾性体の＋ｘ軸断面は、２点鎖線で示すように一様な変形（変位）を行っている。しかし、定数弾性係数の要素で構成した離散要素近似力学系は、直線で示すように質量点によって異なった変位量をもたらしており、連続弾性体の変位と適合性が取れていない。つまり、誤差が生じている。
一方、図７Ｂにおいては、２点鎖線と直線が重なっており、（離散モデルでの単位キューブの変形と連続弾性体の変形が等しい）、離散モデルでの単位キューブの変形、つまり＋ｘ軸断面上の全質量点の変位が連続弾性体の変位を近似しており適合性が取れていることが分かる。つまり、この図面に示す大きさで表すと、この発明を適用することにより誤差がほぼ０になっていることが分かる。

図８Ａ、Ｂは図７Ａ、Ｂと同じ連続弾性体に対して、弾性体を離散化するための空間の分割を不均一にした場合の例である。すなわち、それぞれ異なった寸法を持つ４個の単位キューブで、離散要素近似力学系を構成している。図８Ａ、Ｂに示すＬ１、Ｌ２、Ｌ３、Ｌ４をそれぞれ０．０２ｍ、０．０８ｍ、０．０３ｍ、０．０７ｍとし、他の条件は図７Ａ、Ｂと同様とする。

図８Ａは従来の技術における静的変位の結果で、図７Ａの場合より変位誤差が大きく、つまり、図７Ａより直線と２点鎖線との差が大きくなっており、連続弾性体変形と適合性が取れていない。図７Ａと図８Ａの結果から分かるように、離散要素近似力学系の連続弾性体に対する静的変形適合性においては、離散要素の構成密度が均一でない場合、つまり、分割される複数の単位キューブの各辺が等しくない場合、図７Ａの場合と比べて、変形量の誤差がより大きくなる。

図８Ｂは図８Ａの離散要素近似力学系の構造に、図７Ｂの場合と同じく連続体適合弾性要素を適用した場合、つまりこの発明を適用した静的変形の結果である。図８Ｂより、直線と２点鎖線が重なっており、離散要素の構成密度が均一でない場合でも誤差がほぼ０になっていることが分かる。よって、空間の不均一分割により要素密度（単位空間当たりの離散要素）が大きく変化した離散要素近似系に対してもこの発明は有効に作用しており、連続弾性体の変位を近似し、適合性を確保している。

以上の各実施形態の他、本発明である連続弾性体変形シミュレーション方法は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、連続弾性体変形シミュレーション方法において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。

また、連続弾性体変形シミュレーション方法における処理をコンピュータによって実現する場合、連続弾性体変形シミュレーション方法が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、連続弾性体変形シミュレーション方法における処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Ｄｉｇｉｔａｌ Ｖｅｒｓａｔｉｌｅ Ｄｉｓｃ）、ＤＶＤ−ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ Ｄｉｓｃ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＣＤ−Ｒ（Ｒｅｃｏｒｄａｂｌｅ）／ＲＷ（ＲｅＷｒｉｔａｂｌｅ）等を、光磁気記録媒体として、ＭＯ（Ｍａｇｎｅｔｏ−Ｏｐｔｉｃａｌ ｄｉｓｃ）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Ｅｌｅｃｔｒｏｎｉｃａｌｌｙ Ｅｒａｓａｂｌｅ ａｎｄ Ｐｒｏｇｒａｍｍａｂｌｅ−Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）等を用いることができる。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｅｒｖｉｃｅ Ｐｒｏｖｉｄｅｒ）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、連続弾性体変形シミュレーション方法を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

この発明のハードウェアの構成例を示すブロック図。 連続弾性体における単位キューブを示す図。 離散要素近似力学系における単位キューブを示す図。 離散要素Ａ〜Ｆとノード１〜４による構成されている２次元離散要素近似系を示す図。 図５Ａは剛性マトリクスＫを簡略化した図であり、図５Ｂは図３記載の質量点１、２、３、４を構成する断面に荷重を加えた場合の縮小される剛性マトリクスであり、図５Ｃは図３記載の質量点２、３、５、８を構成する断面に荷重を加えた場合の縮小される剛性マトリクスである。 図６Ａは連続体適合単位キューブＡ、Ｂ、Ｃ、Ｄを表す図であり、図６Ｂは連続体適合単位キューブＡ、Ｂ、Ｃ、Ｄが組み合わされた連続弾性体を示す図である。 連続弾性体が単位キューブに均一に分割された場合、一様断面応力荷重により連続弾性体と離散要素近似力学系の変形を示し、図７Ａは従来の技術を適用した場合の変形を示す図であり、図７Ｂはこの発明を適用した場合の変形を示す図である。 連続弾性体が単位キューブに不均一に分割された場合、一様断面応力荷重により連続弾性体と離散要素近似力学系の変形を示し、図８Ａは従来の技術を適用した場合の変形を示す図であり、図８Ｂはこの発明を適用した場合の変形を示す図である。

Claims (4)

1. 連続弾性体を複数の離散質量点に網状の離散要素を張った離散要素近似力学系モデルの連続弾性体変形シミュレーション方法において、
   連続体面変位計算手段が、隣接している８つの上記離散質量点をそれぞれ頂点とする６面体を単位キューブとし、各単位キューブに対し、連続体としての直交する３軸の各軸方向の荷重に対する連続体面変位を計算する過程と、
   頂点変位計算手段が、各単位キューブに対し、上記離散要素近似力学系における各軸方向の荷重ベクトルに対する各頂点の変位を剛性マトリクスを用いて計算する過程と、上記荷重ベクトルは上記連続体面変位計算過程における対応する単位キューブの同一軸方向の荷重が、各質量点に分割されて印加される構成とされ、
   近似面変位計算手段が、各単位キューブごとに、同一面内の上記頂点変位を平均して近似面変位を求める過程と、
   調整手段が、各単位キューブごとに対応軸面の上記連続体面変位と上記近似面変位との差の絶対値が予め決められた誤差以下になるように上記剛性マトリクスの要素を調整する過程と、を有することを特徴とする連続弾性体変形シミュレーション方法。
2. 請求項１記載の連続弾性体変形シミュレーション方法において、
   共有要素平均手段が、上記単位キューブの共有される離散要素に対し、これらの単位キューブごとの上記調整手段により調整された弾性係数を平均して、上記離散要素の弾性係数として、上記剛性マトリクスの対応する構成要素を修正する過程と、
   修正値頂点変位計算手段が、上記各単位キューブに対し、上記離散要素近似力学系における各軸方向の荷重ベクトルに対する各頂点の変位を上記修正剛性マトリクスを用いて計算する過程と、
   修正面変位計算手段が、各単位キューブごとに、上記修正剛性マトリクスを用いて計算された同一面内の上記頂点変位を平均して修正面変位を求める過程とを、有することを特徴とする連続弾性体変形シミュレーション方法。
3. 請求項１または２の何れかに記載した連続弾性体変形シミュレーション方法の各過程をコンピュータに実行させるための連続弾性体変形シミュレーションプログラム。
4. 請求項３記載の連続弾性体変形シミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体。

Images