(57)【要約】【課題】本発明は、非線形の有限要素法による計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にする構造解析方法と、プログラム、記録媒体、シミュレーション装置を提供することを目的とする。【解決手段】本発明は、非線形の有限要素法による構造解析方法、プログラム、記録媒体、シミュレーション装置であって、３次元要素を特定するプリミティブとして、各３次元要素において他の３次元要素で既に指定されている場合を除いて辺εの組み合わせを指定する辺対Ｐを設定し、辺εの辺ベクトルｐ〜（ε）に対する重み関数ｇ（ε）を計算して、すべての辺対Ｐに関する前記重み関数ｇ（ε）の計算を一括して行うことを主な特徴とする。【選択図】図１
2009301149.tif 000002
【特許請求の範囲】
【請求項１】
複数の３次元要素で分割された対象物モデルの弾性変形を解析するとき、前記対象物モデルの節点νに力ｆが作用した場合に前記対象物モデルがサンブナン・キルヒホフ材料則で変形するとの前提に立って、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法による構造解析方法であって、
前記３次元要素を特定するプリミティブとして、各３次元要素において他の３次元要素で既に指定されている場合を除いて辺εの組み合わせを指定する辺対Ｐを設定し、
前記辺εの辺ベクトルｐ〜（ε）に対する重み関数ｇ（ε）を計算して、すべての辺対Ｐに関する前記重み関数ｇ（ε）の計算を一括して行い、各節点νの位置ｐ（ν）を基に前記辺ベクトルｐ〜（ε）を決定し、
各辺εに関して、前記辺ベクトルｐ〜（ε）と前記重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して前記積の総和をとり、対象物モデルの変形と力ｆの関係を求めることを特徴とする非線形の有限要素法による構造解析方法。
ここに、ｐ〜（ε）は
2009301149.tif 000090
（外字１）を示す。
【請求項２】
複数の四面体要素で分割された対象物モデルの弾性変形を解析するとき、前記対象物モデルの節点νに力ｆが作用した場合に前記対象物モデルがサンブナン・キルヒホフ材料則で変形するとの前提に立って、前記対象物モデルの辺εから任意の辺εａ，辺εｂを取り出し、歪みエネルギーＷτが（数３）に示す関数ζ（ε）を用いた関数ζ（εａ）及び関数ζ（εｂ）とこのζ（εａ），ζ（εｂ）により分解される関数Ｌ（εａ，εｂ）とに基づいて（数１）によって与えられるとして、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法による構造解析方法であって、
各四面体要素を特定するプリミティブとして、前記辺εａ，辺εｂの組み合わせを指定する辺対Ｐを他の四面体要素において既に指定されている場合を除いて設定し、
前記関数Ｌ（εａ，εｂ）を基にして（数２）によって定義される辺εと関係付けられた関数Ｌ（ε）及び辺対Ｐと関係付けられた関数Ｌ（Ｐ）を計算し、
この関数Ｌ（ε）と前記関数ζ（ε）とによってＬ（ε）ζ（ε）を計算すると共に、前記辺対Ｐで指定された辺εP（Ｐ，０），εP（Ｐ，１）ごとに前記関数Ｌ（Ｐ）と前記関数ζ（εP（Ｐ，０）），ζ（εP（Ｐ，１））の各積を計算して前記辺対Ｐに関して総和をとることにより、（数４）で定義される重み関数ｇ（ε）を計算し、
前記対象物モデルにおける前記節点νの空間の中での前記節点νの位置ｐ（ν）を設定して、前記辺εの辺ベクトルｐ〜（ε）をｐ〜（ε）＝ｐ（νｅ（ε，０））−ｐ（νｅ（ε，１））として決定し、
各辺εに関して、前記辺ベクトルｐ〜（ε）と前記重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して前記積の総和をとり、前記節点νに作用する力ｆ（ν）を（数５）に従って計算することを特徴とする非線形の有限要素法による構造解析方法。
【数１】
2009301149.tif 000091
【数２】
2009301149.tif 000092
【数３】
2009301149.tif 000093
【数４】
2009301149.tif 000094
【数５】
2009301149.tif 000095
ここにεＰ（Ｐ，０）は辺対Ｐで指定される一対の辺εの一方（０番目）の辺、εＰ（Ｐ，１）はこの一対の辺εの他方（１番目）の辺であり、νｅ（ε，０）は辺εの一方（０番目）の節点、νｅ（ε，１）は辺εの他方（１番目）の節点、ｐ〜ini（ε）はｐ〜（ε）の初期値を示す。
【請求項３】
請求項２記載の非線形の有限要素法による構造解析方法において、
前記関数Ｌ（Ｐ）と前記辺εＰ（Ｐ，０）、εＰ（Ｐ，１）の辺ベクトルｐ〜（εＰ（Ｐ，０））、ｐ〜（εＰ（Ｐ，０））Ｔから（数６）で決まる関数ＫＰ（Ｐ）を求め、
前記辺εＰ（Ｐ，０），εＰ（Ｐ，１）の各節点に対して（数８）で定まる条件ｃｏｎｄＡ（ε，Ｐ，ｉ，ｊ），ｃｏｎｄＢ（ε，Ｐ，ｉ，ｊ）の真偽値を評価して、これらの条件と前記重み関数ｇ（ε）、前記辺ベクトルｐ〜（ε）、前記関数Ｌ（ε）を使って（数７）により前記対象物モデル全体の接線剛性行列の部分行列Ｋ（ε）を計算することを特徴とする非線形の有限要素法による構造解析方法。
【数６】
2009301149.tif 000096
【数７】
2009301149.tif 000097
【数８】
2009301149.tif 000098
ここに、ｉ，ｊはｉ＝｛０，１｝、ｊ＝｛０，１｝の何れかの整数である。
【請求項４】
コンピュータに、
対象物モデルを分割する３次元要素の節点νに力ｆが作用したとき、前記対象物モデルがサンブナン・キルヒホフ材料則で変形することを前提にして非線形の有限要素法に基づいた計算をさせるために、前記３次元要素の辺εを基に各３次元要素を特定するプリミティブとして、他の３次元要素において既に定義されている組み合わせを除いて前記辺εの一対の組み合わせの関係付けを行う辺対Ｐを設定する手順と、
前記辺εの辺ベクトルｐ〜（ε）に対する重み関数ｇ（ε）を計算して、すべての辺対Ｐに関する前記重み関数ｇ（ε）の計算を一括して行う手順と、
各節点νの位置ｐ（ν）を基に前記辺ベクトルｐ〜（ε）を決定する手順と、
各辺εに関して、前記辺ベクトルｐ〜（ε）と前記重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して前記積の総和をとる手順と、を実行させるためのプログラム。
【請求項５】
コンピュータに、
対象物モデルを分割する四面体要素の節点νに力ｆ（ν）が作用したとき、前記対象物モデルがサンブナン・キルヒホフ材料則で変形することを前提にして非線形の有限要素法に基づいた計算をさせるために、前記対象物モデルの辺εの中から２つの辺εａ，辺εｂを取り出し、歪みエネルギーＷτが（数３）に示す関数ζ（ε）を用いた関数ζ（εａ），関数ζ（εｂ）とこのζ（εａ），ζ（εｂ）により分解される関数Ｌ（εａ，εｂ）とから（数１）によって与えられるとして、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法に基づいて計算させるために、前記四面体要素を特定するプリミティブの１つとして、他の四面体要素において既に定義されている組み合わせを除いて辺εの一対の組み合わせの関係付けを行う辺対Ｐを設定する手順と、
前記辺εを基に（数２）の関数Ｌ（ε）、（数３）の関数ζ（ε）とから重み関数ｇ（ε）を計算する手順と、
辺対Ｐで指定された辺εＰのうち、一方のεＰ（Ｐ，１）について他方のεＰ（Ｐ，０）を使って（数４）の関数Ｌ（Ｐ）とζ（εＰ（Ｐ，０））の積を計算し、辺対Ｐに関してこの積の計算を繰り返し、順次和をとると共に、前記辺εＰの他方のεＰ（Ｐ，０）についてももう一方のεＰ（Ｐ，１）を使って（数４）の関数Ｌ（Ｐ）とζ（εＰ（Ｐ，１））の積を計算し、辺対Ｐに関してこの積Ｌ（Ｐ）ζ（εＰ（Ｐ，１））の計算を繰り返し、順次和をとって重み関数ｇ（ε）を更新する手順と、
辺ベクトルｐ〜（ε）とこの重み関数ｇ（ε）の積を力の初期値とし、各辺εで共通する属性をもつ第１の節点νｅ（ε，０）におけるｆ（νｅ（ε，０））を計算して前記初期値に順に加算すると共に、第２の節点νｅ（ε，１）におけるｆ（νｅ（ε，１））を計算して前記初期値から順に減算し、各節点νにおけるｆ（ν）を（数５）に従って計算する手順と、を実行させるためのプログラム。
【請求項６】
請求項５に記載されたプログラムにおいて、
（数９）でＫ（ε）を計算する手順と、
（数６）でＫｐ（Ｐ）を計算する手順と、
（数１０）でｔ（ν，ε）を定義したとき（数１１）のフラグｓ（Ｐ，ｉ，ｊ）がｓ（Ｐ，ｉ，ｊ）＝１の場合、（数１２）を計算してＫ（ε）に順に加算すると共に、ｓ（Ｐ，ｉ，ｊ）＝−１の場合、（数１３）を計算してＫ（ε）に順に加算する手順と、を実行させるプログラム。
【数９】
2009301149.tif 000099
【数１０】
2009301149.tif 000100
【数１１】
2009301149.tif 000101
【数１２】
2009301149.tif 000102
【数１３】
2009301149.tif 000103
【請求項７】
対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法による構造解析を行うプログラムを格納したコンピュータ読取可能な記憶媒体であって、請求項４乃至請求項６の何れか１項に記載のプログラムを格納すると共に、プリミティブとして少なくとも節点の位置ｐ（ν）、辺、辺対Ｐに関する情報を格納したことを特徴とするコンピュータが読取可能な記憶媒体。
【請求項８】
対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法に基づいた処理を行うシミュレーション装置であって、請求項７記載の記憶媒体で構成された記憶手段と、前記対象物モデルに加える変位量を入力できる入力手段と、前記プログラムを読み出して実行する制御手段とを備えたことを特徴とするシミュレーション装置。
【請求項９】
前記入力手段から前記対象物モデルに変位を加えたとき、前記対象物モデルが変形した状態を表示する表示手段を備えたことを特徴とする請求項８記載のシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象物モデルを近似的に表現する非線形の有限要素法による構造解析方法であって、とくに大きな変形に対しても精度高く対象物モデルを表現でき、動的な計算においても３次元モデルの実時間の動作を可能にする構造解析方法と、プログラムと、そのプログラムを格納した記録媒体、さらに対象物モデルを近似的に表現し、大きな変形に対しても高精度で表現できるシミュレーション装置に関する。
【背景技術】
【０００２】
近年、外科手術の技術向上などを目的として、手術用のシミュレーション装置が開発されている。このシミュレーション装置で使用される臓器モデルとしては、有限要素法を用いたモデルが利用されることが多い。有限要素法は、モデルを細かな多面体要素に分割し、要素の節点に働く力と節点の変位の関係を簡単な近似式に置き換えることで、モデル全体の変形を解析する方法である。
【０００３】
ところで、有限要素法としては、要素内の変位が微小という前提をおいて、力と変位の関係を線形化したモデルが多用される。この線形の有限要素法では線形のコーシー歪みを用いて要素内の歪みを表すことになる。歪みの２次の項は微小という仮定を必要とする。しかし、これでは回転などの大きな変形に対応できないため、大きな変形を表現できないという課題があった。
【０００４】
そこで、従来、線形の有限要素法を大きな変形にも対応できるようにする提案がなされている。１つの方法は、線形性を保てる範囲の変形を行い、逐次接線剛性行列を更新し、大きな変形を微小変形の重ね合わせとして扱うものである。また、高周波成分を取り除くことで、非線形の変形モデルを線形の変形モデルに置き換える手法も提案されている。
【０００５】
もう１つの方法は非線形の有限要素法を用いる方法である。すなわち、コーシー歪みに代わって、２次の項を考慮したグリーン・ラグランジュ歪みを用いる（例えば非特許文献１）。これにより、回転や変位などの大きな変形にも対応できるようになる。しかし、このグリーン・ラグランジュ歪みを使ったモデルは、変形を精度よく表現することはできるが、歪みの２次の項が含まれるため方程式が複雑化し、節点に働く力と変位との関係を表す連立方程式を簡単に解くことが困難になる。
【０００６】
この解法として、ニュートン・ラプソン法などの反復計算を繰り返す方法もあるが、１ステップあたりの計算時間がかかりすぎる。すなわち、計算量が増加することで、実時間で対応することが困難になる。この計算量の低減が非線形の有限要素法における最大の問題点である。
【０００７】
そこで、事前計算によって求めた変形モードの重ね合わせを用いることで、実計算の計算時間の短縮が図られている。例えば、モデルの各点の変位をエネルギーに結びつける力学テンソルを多項式に展開し、大きな変位に対して、並進移動による変形、回転運動で均一な変形、均一ではないが比較的広い範囲の歪みを有する変形、局所的な大きな変形、の変形モードを設けて対処する提案がある（特許文献１参照）。また、モデルの「骨格」を定義した上で、その骨格の周囲の要素を局所的な線形モデルとして近似することで計算量を減らすことも試みられている。
【０００８】
さて、有限要素法の計算には逆演算なしに演算する陽解法と、逆演算を含む形で演算する陰解法、の２通りの計算方法がある。陰解法はモデルの挙動のステップ時間の長さに対する影響が小さく、解が発散しにくい。これに対し陽解法は発散する可能性はあるが計算量が少なくて済む。なお、陰解法には静的と動的の２つの場合がある。前者はモデルが運動を考慮に入れない場合であり、後者は運動を考慮に入れる場合である。陽解法は通常はモデルの運動を考慮に入れて行われる。
【０００９】
動的陰解法の１つ目は、特許文献１のように事前計算により変形モードを求めておき、それを重ね合わせることができると考える手法であり、モード単位にすることで非線形歪みの計算量を削減する。２つ目は非線形歪みを適用する範囲を限定する手法であり、モデルの「骨格」を定義し、その骨格の周囲では要素を線形モデルと近似することで、非線形として計算する範囲を限定して計算量を減らす。
【００１０】
また、動的陽解法を用いた非線形の有限要素法の解法も提案されている（非特許文献２，３）。この動的陽解法は、ステップ時間の長さに影響されやすいと言う欠点もあるが、事前計算や反復計算が必要でなく、演算の並列化が容易という利点もある。非特許文献２の動的陽解法は、グリーン・ラグランジュ歪みに基づいており、節点に働く力を節点の変位の３次多項式として表している。この方法は線形の有限要素法に比べて、数倍の計算時間がかかる。
【００１１】
【特許文献１】特開２００３−３１６８３３号公報
【非特許文献１】H.Delingette, 「Triangular Springs for Modeling Nonlinear Membranecs,IEEE Transactions on visualization and computer graphics,vol14,No.2,2008,p.329-p.341
【非特許文献２】G.Picinbomo, H.Delingette, N.Ayache, 「Non-linear anisotoropicelasticity for real-time surgery simulation」, Graphical Models,65(5),2003, p.305-p.321
【非特許文献３】田淵裕章、菊植亮、山本元司、「グリーン・ラグランジュひずみに基づいた実時間有限要素法の省計算量表現」、第１２回日本バーチャルリアリティ学会講演論文集、講演番号１Ｄ４−１、２００７年９月１９日
【発明の開示】
【発明が解決しようとする課題】
【００１２】
以上説明したように、従来の線形の有限要素法においては、２次の項を省略した線形のコーシー歪みを用いて要素内の歪みを表している。これに対して非線形の有限要素法では、２次の項を考慮したグリーン・ラグランジュ歪みを用いる。これにより、回転などの大きな変形にも対応できるようになる。
【００１３】
このグリーン・ラグランジュ歪みを使った非線形の有限要素法モデルは、モデルの変形を高精度で表現できるが、力を変位との関係を表す連立方程式を簡単には解くことができない。２次の項が含まれるため方程式が複雑化するためである。また、計算量が増加することで、実時間での対応が困難になる。この計算量の低減が非線形の有限要素法における課題であった。
【００１４】
非特許文献１の手法は、最大２次元の解析であって、非線形のグリーン・ラグランジュ歪みに基づき、三角形要素を用い、節点に働く力を節点の変位の３次の多項式として表し、これを三角形要素の辺間の内角を使うことで計算している。しかし、２次元の解析では３次元と比べて利用価値が低く、利用価値を高めるためには２次元の解析を３次元に拡張しなければならない。しかし、この解析方法は三角形要素の内角を使うため拡張性に問題があるものであった。
【００１５】
本発明者らは、非特許文献３において、要素の辺を使用することで、非特許文献２の有限要素法の式を別の式に導き直し、計算量を減らせることを開示した。有限要素法のプリミティブには、節点、辺、面、要素の形状が考えられるが、節点と辺を使っての計算量の低減であった。
【００１６】
非特許文献２の手法もグリーン・ラグランジュ歪みに基づいており、節点に働く力を変位の３次多項式として表している。ただ、コーシー歪みに基づいた線形の有限要素法に比べて、５倍程度の計算時間がかかる。動的陽解法を使った場合、実時間での動作は困難になる。そこで、非特許文献２においてはグリーン・ラグランジュ歪みを用いて計算する範囲を限定し、大部分をコーシー歪みで計算する計算手法を使用している。
【００１７】
このように従来の非線形の有限要素法は仮定や条件を設けて近似したり、複数の手法を寄せ集めて計算量を低減しているが、本格的、飛躍的な計算量の低減にはなっておらず、より高速の演算を行うためには新たな視点による本格的な計算量の低減が望まれる。
【００１８】
そして、手術用のシミュレーション装置だけでなく、３Ｄモデル教育やコンピュータゲームなど、インタラクティブなアプリケーションにとって、実時間での動作は不可欠である。モデルの動作や形状が不自然になったり、歪みができたりしてはこうした装置、あるいはアニメーションにも利用できない。
【００１９】
そこで、本発明は、非線形の有限要素法による計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にする構造解析方法と、非線形の有限要素法による計算量を大幅に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にするプログラムと、そのプログラムを格納した記録媒体、さらに、非線形の有限要素法による計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にするシミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【００２０】
本発明の非線形の有限要素法による構造解析方法は、複数の３次元要素で分割された対象物モデルの弾性変形を解析するとき、前記対象物モデルの節点νに力ｆが作用した場合に前記対象物モデルがサンブナン・キルヒホフ材料則で変形するとの前提に立って、対象物モデル全体の変形と力ｆとの関係を解析する構造解析方法であって、前記３次元要素を特定するプリミティブとして、各３次元要素において他の３次元要素で既に指定されている場合を除いて辺εの組み合わせを指定する辺対Ｐを設定し、前記辺εの辺ベクトルｐ〜（ε）に対する重み関数ｇ（ε）を計算して、すべての辺対Ｐに関する前記重み関数ｇ（ε）の計算を一括して行い、各節点νの位置ｐ（ν）を基に前記辺ベクトルｐ〜（ε）を決定し、各辺εに関して、前記辺ベクトルｐ〜（ε）と前記重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して前記積の総和をとり、対象物モデルの変形と力ｆの関係を求めることを主な特徴とする。ここに、ｐ〜（ε）は
2009301149.tif 000003
（外字１）を示す。
【発明の効果】
【００２１】
本発明の非線形の有限要素法による構造解析方法、プログラム、記録媒体、シミュレーション装置によれば、サンブナン・キルヒホフ材料則に基づいた非線形の有限要素法によって構造解析するときに計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【発明を実施するための最良の形態】
【００２２】
上記課題を解決するための本発明の第１の形態は、複数の３次元要素で分割された対象物モデルの弾性変形を解析するとき、前記対象物モデルの節点νに力ｆが作用した場合に前記対象物モデルがサンブナン・キルヒホフ材料則で変形するとの前提に立って、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法による構造解析方法であって、前記３次元要素を特定するプリミティブとして、各３次元要素において他の３次元要素で既に指定されている場合を除いて辺εの組み合わせを指定する辺対Ｐを設定し、前記辺εの辺ベクトルｐ〜（ε）に対する重み関数ｇ（ε）を計算して、すべての辺対Ｐに関する前記重み関数ｇ（ε）の計算を一括して行い、各節点νの位置ｐ（ν）を基に前記辺ベクトルｐ〜（ε）を決定し、各辺εに関して、前記辺ベクトルｐ〜（ε）と前記重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して前記積の総和をとり、対象物モデルの変形と力ｆの関係を求めることを特徴とする非線形の有限要素法による構造解析方法である。
ここに、ｐ〜（ε）は
2009301149.tif 000004
（外字２）を示す。
【００２３】
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によって構造解析するときに、辺対Ｐをプリミティブとして採用することにより各要素間において重複する計算を省くことができ、計算量を低減することが可能になり、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００２４】
本発明の第２の形態は、複数の四面体要素で分割された対象物モデルの弾性変形を解析するとき、前記対象物モデルの節点νに力ｆが作用した場合に前記対象物モデルがサンブナン・キルヒホフ材料則で変形するとの前提に立って、前記対象物モデルの辺εから任意の辺εａ，辺εｂを取り出し、歪みエネルギーＷτが（数１６）に示す関数ζ（ε）を用いた関数ζ（εａ）及び関数ζ（εｂ）とこのζ（εａ），ζ（εｂ）により分解される関数Ｌ（εａ，εｂ）とに基づいて（数１４）によって与えられるとして、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法による構造解析方法であって、各四面体要素を特定するプリミティブとして、前記辺εａ，辺εｂの組み合わせを指定する辺対Ｐを他の四面体要素において既に指定されている場合を除いて設定し、前記関数Ｌ（εａ，εｂ）を基にして（数１５）によって定義される辺εと関係付けられた関数Ｌ（ε）及び辺対Ｐと関係付けられた関数Ｌ（Ｐ）を計算し、この関数Ｌ（ε）と前記関数ζ（ε）とによってＬ（ε）ζ（ε）を計算すると共に、前記辺対Ｐで指定された辺εP（Ｐ，０），εP（Ｐ，１）ごとに前記関数Ｌ（Ｐ）と前記関数ζ（εP（Ｐ，０）），ζ（εP（Ｐ，１））の各積を計算して前記辺対Ｐに関して総和をとることにより、（数１７）で定義される重み関数ｇ（ε）を計算し、前記対象物モデルにおける前記節点νの空間の中での前記節点νの位置ｐ（ν）を設定して、前記辺εの辺ベクトルｐ〜（ε）をｐ〜（ε）＝ｐ（νｅ（ε，０））−ｐ（νｅ（ε，１））として決定し、各辺εに関して、前記辺ベクトルｐ〜（ε）と前記重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して前記積の総和をとり、前記節点νに作用する力ｆ（ν）を（数１８）に従って計算することを特徴とする非線形の有限要素法による構造解析方法である。
【００２５】
【数１４】
2009301149.tif 000005
【００２６】
【数１５】
2009301149.tif 000006
【００２７】
【数１６】
2009301149.tif 000007
【００２８】
【数１７】
2009301149.tif 000008
【００２９】
【数１８】
2009301149.tif 000009
ここにεＰ（Ｐ，０）は辺対Ｐで指定される一対の辺εの一方（０番目）の辺、εＰ（Ｐ，１）はこの一対の辺εの他方（１番目）の辺であり、νｅ（ε，０）は辺εの一方（０番目）の節点、νｅ（ε，１）は辺εの他方（１番目）の節点、ｐ〜ini（ε）はｐ〜（ε）の初期値を示す。
【００３０】
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によって構造解析するときに、辺対Ｐをプリミティブとして採用することにより各要素間において重複する計算を省くことができ、対象物モデルの各節点に働く力を計算するときの計算量を大幅に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００３１】
本発明の第３の形態は、第２の形態に従属する形態であって、前記関数Ｌ（Ｐ）と前記辺εＰ（Ｐ，０）、εＰ（Ｐ，１）の辺ベクトルｐ〜（εＰ（Ｐ，０））、ｐ〜（εＰ（Ｐ，０））Ｔから（数１９）で決まる関数ＫＰ（Ｐ）を求め、前記辺εＰ（Ｐ，０），εＰ（Ｐ，１）の各節点に対して（数２１）で定まる条件ｃｏｎｄＡ（ε，Ｐ，ｉ，ｊ），ｃｏｎｄＢ（ε，Ｐ，ｉ，ｊ）の真偽値を評価して、これらの条件と前記重み関数ｇ（ε）、前記辺ベクトルｐ〜（ε）、前記関数Ｌ（ε）を使って（数２０）により前記対象物モデル全体の接線剛性行列の部分行列Ｋ（ε）を計算することを特徴とする非線形の有限要素法による構造解析方法である。
【００３２】
【数１９】
2009301149.tif 000010
【００３３】
【数２０】
2009301149.tif 000011
【００３４】
【数２１】
2009301149.tif 000012
ここに、ｉ，ｊはｉ＝｛０，１｝、ｊ＝｛０，１｝の何れかの整数である。
【００３５】
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によって構造解析するときに、辺対Ｐをプリミティブとして採用することにより各要素間で重複する計算を省くことができ、対象物モデルの剛性を計算するときの計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００３６】
本発明の第４の形態は、コンピュータに、対象物モデルを分割する３次元要素の節点νに力ｆが作用したとき、前記対象物モデルがサンブナン・キルヒホフ材料則で変形することを前提にして非線形の有限要素法に基づいた計算をさせるために、前記３次元要素の辺εを基に各３次元要素を特定するプリミティブとして、他の３次元要素において既に定義されている組み合わせを除いて前記辺εの一対の組み合わせの関係付けを行う辺対Ｐを設定する手順と、前記辺εの辺ベクトルｐ〜（ε）に対する重み関数ｇ（ε）を計算して、すべての辺対Ｐに関する前記重み関数ｇ（ε）の計算を一括して行う手順と、各節点νの位置ｐ（ν）を基に前記辺ベクトルｐ〜（ε）を決定する手順と、各辺εに関して、前記辺ベクトルｐ〜（ε）と前記重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して前記積の総和をとる手順と、を実行させるためのプログラムである。
【００３７】
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によって構造解析するときに、辺対Ｐをプリミティブとして採用することにより各要素間において重複する計算を省くことができ、対象物モデルの各節点に働く力を計算するときの計算量を大幅に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００３８】
本発明の第５の形態は、コンピュータに、対象物モデルを分割する四面体要素の節点νに力ｆ（ν）が作用したとき、前記対象物モデルがサンブナン・キルヒホフ材料則で変形することを前提にして非線形の有限要素法に基づいた計算をさせるために、前記対象物モデルの辺εの中から２つの辺εａ，辺εｂを取り出し、歪みエネルギーＷτが（数１６）に示す関数ζ（ε）を用いた関数ζ（εａ），関数ζ（εｂ）とこのζ（εａ），ζ（εｂ）により分解される関数Ｌ（εａ，εｂ）とから（数１４）によって与えられるとして、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法に基づいて計算させるために、前記四面体要素を特定するプリミティブの１つとして、他の四面体要素において既に定義されている組み合わせを除いて辺εの一対の組み合わせの関係付けを行う辺対Ｐを設定する手順と、前記辺εを基に（数１５）の関数Ｌ（ε）、（数１８）の関数ζ（ε）とから重み関数ｇ（ε）を計算する手順と、辺対Ｐで指定された辺εＰのうち、一方のεＰ（Ｐ，１）について他方のεＰ（Ｐ，０）を使って（数１７）の関数Ｌ（Ｐ）とζ（εＰ（Ｐ，０））の積を計算し、辺対Ｐに関してこの積の計算を繰り返し、順次和をとると共に、前記辺εＰの他方のεＰ（Ｐ，０）についてももう一方のεＰ（Ｐ，１）を使って（数１７）の関数Ｌ（Ｐ）とζ（εＰ（Ｐ，１））の積を計算し、辺対Ｐに関してこの積Ｌ（Ｐ）ζ（εＰ（Ｐ，１））の計算を繰り返し、順次和をとって重み関数ｇ（ε）を更新する手順と、辺ベクトルｐ〜（ε）とこの重み関数ｇ（ε）の積を力の初期値とし、各辺εで共通する属性をもつ第１の節点νｅ（ε，０）におけるｆ（νｅ（ε，０））を計算して前記初期値に順に加算すると共に、第２の節点νｅ（ε，１）におけるｆ（νｅ（ε，１））を計算して前記初期値から順に減算し、各節点νにおけるｆ（ν）を（数１８）に従って計算する手順と、を実行させるためのプログラムである。
【００３９】
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によって構造解析するときに、辺対Ｐをプリミティブとして採用することにより各要素間で重複する計算を省くことができ、対象物モデルの各節点に働く力を計算するときの計算量を大幅に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００４０】
本発明の第６の形態は、第５の形態に従属する形態であって、（数２２）でＫ（ε）を計算する手順と、（数１９）でＫｐ（Ｐ）を計算する手順と、（数２３）でｔ（ν，ε）を定義したとき（数２４）のフラグｓ（Ｐ，ｉ，ｊ）がｓ（Ｐ，ｉ，ｊ）＝１の場合、（数２５）を計算してＫ（ε）に順に加算すると共に、ｓ（Ｐ，ｉ，ｊ）＝−１の場合、（数２６）を計算してＫ（ε）に順に加算する手順と、を実行させるプログラム
【００４１】
【数２２】
2009301149.tif 000013
【００４２】
【数２３】
2009301149.tif 000014
【００４３】
【数２４】
2009301149.tif 000015
【００４４】
【数２５】
2009301149.tif 000016
【００４５】
【数２６】
2009301149.tif 000017
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によって構造解析するときに、辺対Ｐをプリミティブとして採用することにより各要素間において重複する計算を省くことができ、対象物モデルの剛性を計算するときの計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００４６】
本発明の第７の形態は、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法による構造解析を行うプログラムを格納したコンピュータ読取可能な記憶媒体であって、第４乃至第５の何れか１の形態のプログラムを格納すると共に、プリミティブとして少なくとも節点の位置ｐ（ν）、辺、辺対Ｐに関する情報を格納したことを特徴とするコンピュータが読取可能な記憶媒体である。
【００４７】
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によって構造解析するときに、辺対Ｐをプリミティブとして採用することにより各要素間において重複する計算を省くことができ、計算量を大幅に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００４８】
本発明の第８の形態は、対象物モデル全体の変形と力ｆとの関係を解析する非線形の有限要素法に基づいた処理を行うシミュレーション装置であって、第７の形態の記憶媒体で構成された記憶手段と、前記対象物モデルに加える変位量を入力できる入力手段と、前記プログラムを読み出して実行する制御手段とを備えたことを特徴とするシミュレーション装置である。
【００４９】
この構成によって、対象物モデルがサンブナン・キルヒホフ材料則で変形する非線形の有限要素法によってシミュレーションするときに、辺対Ｐをプリミティブとして採用することにより各要素間で重複する計算を省くことができ、計算量を大幅に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【００５０】
本発明の第９の形態は、第８の形態に従属する形態であって、入力手段から対象物モデルに変位を加えたとき、対象物モデルが変形した状態を表示する表示手段を備えたことを特徴とするシミュレーション装置である。
【００５１】
この構成によって、対象物モデルに変位を加えたとき、表示手段に対象物モデルの変形状態を表示することができる。
【００５２】
（実施の形態１）
以下、非線形の有限要素法による計算量を本格的に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にする本発明の実施の形態１におけるプログラムと、そのプログラムを格納した記録媒体、画像形成装置について説明をする。
【００５３】
図１は本発明の実施の形態１における非線形の有限要素法の構造解析方法が前提とする対象物モデルのメッシュの説明図、図２は四面体要素の一面を共通にして接続した各プリミティブの説明図、図３は本発明の実施の形態１の非線形の有限要素法による構造解析方法で力ｆを計算するときのフローチャート、図４は本発明の実施の形態１におけるプログラムの各節点の力を演算するアルゴリズムの説明図、図５は辺対が共通の節点をもつ辺か否かで２通りのタイプに分かれることの説明図、図６は本発明の実施の形態１における非線形の有限要素法による構造解析方法で対象物モデル全体の接線剛性行列Ｋを計算するときのフローチャート、図７は本発明の実施の形態１におけるプログラムの接線剛性行列の非対角成分を演算するアルゴリズムの説明図、図８は本発明の実施の形態１におけるプログラムの接線剛性行列の対角成分を演算するアルゴリズムの説明図、図９は本発明の実施の形態１におけるプログラムの節点に作用する力と対象物モデル全体の接線剛性行列を演算するアルゴリズムの説明図、図１０は本発明の実施の形態１におけるプログラムの一連のアルゴリズムを演算するアルゴリズムの説明図、図１１は本発明の実施の形態１におけるシミュレーション装置の構成図、図１２は本発明の実施の形態１におけるプログラムの後退（陰的）オイラー法で数値計算を行うアルゴリズムの説明図、図１３（ａ）は板状モデルの初期の状態を示す図、図１３（ｂ）は板状モデルを大きく変形して折り曲げた状態を示す図、図１３（ｃ）は板状モデルに大きなねじりを加えた状態を示す図の説明図である。
【００５４】
［力ｆと辺ベクトル
2009301149.tif 000018
（外字２）の計算］
図１に示すように、Ｍは空間に存在するモデル（本発明の対象物モデル）の中に形成されたメッシュであり、τはメッシュＭによって構成される四面体要素（本発明における３次元要素）、Ｓはこの四面体要素τを構成する面、εはその辺、νは四面体要素τの節点である。
【００５５】
従って、メッシュＭ中の１つの四面体要素τには、節点νが４個存在し、これによって４個の面Ｓが形成されると共に、辺εが６個形成される。辺εは一対の節点ν（例えばν＝ａとν＝ｂ）によって決定されることになる。以下、例えば節点ａを要素とする辺εをεａ、また、節点ｂを要素とする辺εをεｂと表記する。この四面体要素τ、面Ｓ、辺ε、節点νはモデルを構築するプリミティブである。但し、本発明の要点はモデルをメッシュＭでどのように分割するかと言うより、モデルがメッシュＭによって分割された状態でどのように計算を進めれば演算の高速化が図れるかと言う点にある。従って、以下、モデルが所与の３次元要素（四面体要素τ）で分割されているものとして説明する。
【００５６】
さて、ここでメッシュＭの中にＮ個の節点νが存在するものとし、ｐ（ν）を節点νの現在の位置ベクトルとし（ｐ（ν）∈Ｒ３）、ベクトルｆ（ν）を節点νに働く力（ｆ（ν）∈Ｒ３）とすると、Ｎ個のνに対して力ｆはｐの関数としてｆ＝ｆ（ｐ）のように記述できる。ここで、ｆ∈Ｒ３Ｎであり、ｐ∈Ｒ３Ｎである。サンブナン・キルヒホフ材料則が支配する物体の場合は、ｆ（ｐ）はｐの３次多項式で表される。
【００５７】
次に（数２７）によって行列Ｋ（ｐ）を定義する。このＫ（ｐ）は対象物モデル全体の接線剛性行列であり、キルヒホフの法則に従う物体の場合は、位置ベクトルｐの２次関数になる。
【００５８】
【数２７】
2009301149.tif 000019
ところで、節点νの初期の位置ベクトルをｐini（ν）、現在の位置ベクトルをｐ（ν）とすると、グリーン・ラグランジュ歪みは（数２８）のように歪みテンソルＥとして表現できる。歪みテンソルＥはテンソル積ＦＴＦに関係付けられる。ここでＦは変位勾配テンソルを表す行列である。なお、ＦＴはＦの転置行列であり、Ｉｎはn次の単位行列である。そして、サンブナン・キルヒホフ材料則（以下、サンブナン・キルヒホフの法則とも言う）は一般式の（数２９）で定義される。
【００５９】
【数２８】
2009301149.tif 000020
【００６０】
【数２９】
2009301149.tif 000021
簡単のため、（数３０）により歪みベクトルεを定義する。（数３０）のエレメントは（数３１）で記述されたテンソルＥから取り出されたものである。テンソルＥは対称なテンソルである。
【００６１】
【数３０】
2009301149.tif 000022
【００６２】
【数３１】
2009301149.tif 000023
サンブナン・キルヒホフの法則が適用されるような等方均質の場合は、（数３０）の歪みベクトルεを使って歪みエネルギー密度ｗは（数３２）で表される。ここに、εＴは歪みベクトルεの転置ベクトルであり、Ｄは（数３３）で定義される行列である。λ、μはラーメ定数である。これをヤング率Ｅ、ポアソン比νを使って表すと、λ＝Ｅν／（１＋ν）／（１−２ν）、μ＝Ｅ／２／（１＋ν）となる。
【００６３】
【数３２】
2009301149.tif 000024
【００６４】
【数３３】
2009301149.tif 000025
次に、辺εを使った新しい解析による歪みエネルギーの誘導について説明する。これは、非線形のグリーン・ラグランジュ歪みを前提にしたときでも、メッシュＭの中の各四面体要素τの歪みエネルギーが四面体要素τの１組の辺εに対して分解されることであり、本発明の実施の形態１における非線形の有限要素法はこの性質を利用して隣接する四面体要素τ間で重畳する演算をすることなく計算を進めることができる。なお、従来はメッシュＭの所定形状の要素単位で演算をし、要素毎（エレメント・バイ・エレメント）へと演算を行っていた。
【００６５】
まず、要素内でサンブナン・キルヒホフの法則に従う物体中の四面体要素τを考える。この四面体要素τには６個の辺εと４個の節点νが存在する。ｐini（ν）∈Ｒ３が節点νの初期の位置とすると、辺の長さ
2009301149.tif 000026
（外字２）、
2009301149.tif 000027
（外字３）は（数３４）、（数３５）で定義される。ここで、νｅ（ε，ｊ）は辺εのｊ番目の節点であり、辺εは一対の節点からなるためｊ＝０，１である。
【００６６】
【数３４】
2009301149.tif 000028
【００６７】
【数３５】
2009301149.tif 000029
さて、辺ベクトル勾配行列Ｆが四面体要素τの中で定数行列Ｆτの場合、（数２９）に従って、
2009301149.tif 000030
（外字２）と
2009301149.tif 000031
（外字３）は（数３６）のような関係をもつ。
【００６８】
【数３６】
2009301149.tif 000032
（数２８）から四面体要素τの中のＥが変位勾配行列Ｆにより（数３７）によって得られるが、計算を簡単にするため
2009301149.tif 000033
（外字４）を（数３８）によって定義すると、（数３９）の関係がある。
【００６９】
【数３７】
2009301149.tif 000034
【００７０】
【数３８】
2009301149.tif 000035
【００７１】
【数３９】
2009301149.tif 000036
四面体要素τの６個の辺をεｉ（ｉ＝０，１，２，３，４，５）とし、Ｑτ、ζτ、ζ（ε）、Ｈ（ε）を（数４０）、（数４１）、（数４２）、（数４３）で定義すると、Ｑτετ＝ζτ／２となる。そして、Ｑτには逆行列Ｑτ−１が存在するから、ετ＝（１／２）Ｑτ−１ζτが得られる。このετと（数３２）とから四面体要素τの歪みエネルギーＷτを表す（数４４）が導かれる。ここに、ζ（ε）は辺εの特徴量であり、これはメッシュを構成する四面体要素τの形状に依存しない量である。これに対しＬτは（数４５）で定義され、メッシュを構成する四面体要素τに依存する特徴量である。また、Ｃτｉｎｉは四面体要素τの初期有していた体積の６倍の値を示す。
【００７２】
【数４０】
2009301149.tif 000037
【００７３】
【数４１】
2009301149.tif 000038
【００７４】
【数４２】
2009301149.tif 000039
【００７５】
【数４３】
2009301149.tif 000040
【００７６】
【数４４】
2009301149.tif 000041
【００７７】
【数４５】
2009301149.tif 000042
さらに、（数４０）、（数３３）（数４５）を考慮してＬτの（ｉ，ｊ）番目の要素をＬτ（εｉ，εｊ）とすると、（数４４）は（数４６）のように書き直せる。これは四面体要素τの歪みエネルギーであり、本発明の非線形の有限要素法を計算する上で基本となる式である。そして、メッシュＭ全体の歪みエネルギーＷはメッシュＭ全体で四面体要素τに対する歪みエネルギーＷτの総和をとったものであるから、Ｌ（εａ，εｂ）＝Στ∈MＬτ（εａ，εｂ）であり、（数４７）のように記述できることになる。
【００７８】
【数４６】
2009301149.tif 000043
【００７９】
【数４７】
2009301149.tif 000044
この（数４６）の意義を説明すると、非線形のグリーン・ラグランジュ歪みを前提にしたエネルギーＷτが、一対の辺εａ、εｂの各関数であるζ（εａ）、ζ（εｂ）、さらに辺εａ，εｂの関数であるＬ（εａ，εｂ）によって、きわめてきれいな形で３つの関数の積として分解された、すなわち簡単に言えば、いわば因数分解できたことになる。この行列Ｌ（εａ，εｂ）は、εａ＝εｂとなる対角成分と、εａ≠εｂとなる対角成分以外の要素を含んでいる。対角成分はεａ＝εｂであり、これは辺自身によって起こる辺への影響を示し、非対角成分はεａ≠εｂの辺εｂが辺εａに与える影響を示すものである。ここで、力ｆ（ν）は歪みエネルギーＷをｐ（ν）Ｔで偏微分することにより得られ、さらに辺ベクトル
2009301149.tif 000045
（外字２）に対する重み関数ｇ（ε）を（数４８）、ｔ（ν，ε）を（数４９）と定義すると、力ｆ（ν）は（数５０）で与えられる。なお、ｔ（ν，ε）は各辺εにおいて共通する属性（例えば所定方向にみたときの大小）をもつ節点νに対して共通して正負を付与するものである。それ以外の節点では０となる。
【００８０】
【数４８】
2009301149.tif 000046
【００８１】
【数４９】
2009301149.tif 000047
【００８２】
【数５０】
2009301149.tif 000048
（数４０）、（数４５）、Ｌ（εａ，εｂ）＝Στ∈MＬτ（εａ，εｂ）から、Ｌ（εａ，εｂ）は辺εａ，εｂが共通の四面体要素τに属しているときのみ「０」でない値をもつ。従って、共通する四面体要素τに属している辺εのうち、一対の辺εａ，εｂと関係付けられた幾何学的特徴変数、すなわち新しいプリミティブを定義すれば、エレメント・バイ・エレメントで計算を進めるのではなく、この新しいプリミティブを使ってエレメント間で情報を共有しながら、プリミティブごとに計算を進めることが可能になり、重畳する演算は利用し合うことによって計算量を減らすことができる。
【００８３】
ここで、この四面体要素τの６個ある辺から構成される２７個の辺の組み合わせにおいて、一対の辺εａ，εｂの組み合わせをＴＳＥＰ（tetrahedron-sharing edge pairs）と称することにし、四面体要素辺対（本発明における辺対、以下、辺対）と言う。すると、辺対Ｐは共通する１個のメッシュ内で定義され、図２に示すように四面体要素τの場合は１５対の組み合わせが存在し、辺対Ｐの数は１５である。四面体要素τのモデルは４個の節点νを有しており、この２つの節点νの組から６個の辺εができ、３個の節点νの組から４個の面Ｓ、また４個の節点νから１個の四面体要素τができるが、この辺εａ，εｂからその組み合わせ（εａ，εｂ）を示す辺対Ｐの数は１５となる。
【００８４】
しかしながら、ある四面体要素τ０に対して別の四面体要素τ１を追加するときには、追加される四面体要素τ１の辺対Ｐの数は必ずしも一定値１５ではなくなる。図２に示すように、四面体要素τ０に１個の面Ｓを共有して別の四面体要素τ１を追加した場合、四面体要素τ１の辺対Ｐの数は１２（１２対）となる。すなわち、四面体要素τ０、四面体要素τ１間で共通する面Ｓの３個の節点、３個の辺ε、３個の辺対Ｐは、既存の四面体要素τ０の要素として既に定義されており、追加する四面体要素τ１では定義する必要がない。四面体要素τ１で新たに定義しなければならないのは１２個（１２対）であり、四面体要素τ０、τ１全体で必要な辺対Ｐの数は２７（２７対）である。
【００８５】
なお、図２の四面体要素τ０、τ１においては、辺対ＰをＰ０〜Ｐ２７まで通し番号を付して記載しており、辺ε３，辺ε４，辺ε５が四面体要素τ０、τ１間で共有される辺であって、四面体要素τ０、τ１間で共有される辺対ＰはＰ１２，Ｐ１３，Ｐ１４である。これは四面体要素τ０、τ１が図２のような１個の面Ｓで接続されていることを一義的に示し、（辺ε３，辺ε４）、（辺ε４，辺ε５）、（辺ε３，辺ε５）に関する演算が四面体要素τ０、τ１間で共有でき、一方で一度計算したら演算結果を共有可能なことを示す。
【００８６】
同様に、既存の四面体要素τ０、τ１に更に別の１個の四面体要素τ２を追加した場合に、四面体要素τ２の追加による辺対Ｐの追加は新たに９（９対）であり（全辺対Ｐは３６）、既存の四面体要素τ０、τ１、τ２に四面体要素τ３を追加した場合に、四面体要素τ２による四面体要素辺対Ｐの追加は６（６対）となる（全四面体要素辺対Ｐは４２）。そして、既存の四面体要素τ０、τ１、τ２、τ３に１個の四面体要素τ４を追加した場合には、辺対Ｐの追加は０（０対）となる（全辺対Ｐは４２）。
【００８７】
そして、実施の形態1の場合、既存の組み合わせ（εａ，εｂ）を使って追加される辺εを計算するので、エレメント・バイ・エレメントの計算では繰り返される計算を共有することができ、計算速度が速くなる。なお、以上説明したのは３次元要素の中でも四面体要素であるが、辺対Ｐを３次元要素のプリミティブとして採用して、隣接する要素間で演算結果を共有すれば高速演算が可能になるのは、他の３次元要素でも同様である。
【００８８】
しかも、この新たなプリミティブの辺対Ｐを使えば、Ｌ（εａ，εｂ）の対角成分に関する辺εの計算と、非対角成分に関する辺対Ｐ関係の計算に分けて計算することができる。以下その説明を行う。
【００８９】
辺対Ｐのｉ番目（ｉ＝０，１）の辺εをεＰ（Ｐ，ｉ）と表記することにする。例えば、辺対Ｐが（εａ，εｂ）の組み合わせの場合、εａがεＰ（Ｐ，０）と表現でき、εｂがεＰ（Ｐ，１）と表現できる。Ｌ（εａ，εｂ）が「０」でない限り、これは対角成分（εａ＝εｂ）を示すＬ（ε、ε）と、非対角成分（εａ≠εｂ）を示すＬ（Ｐ）とからなる。これを（数５１）のよう表記する。
【００９０】
【数５１】
2009301149.tif 000049
これを使って（数４８）の重み関数ｇ（ε）を書き換えると、（数５２）のようになる。（数２６）の第１項は各辺εの対角成分による重み成分で、第２項は辺対Ｐのうち一方（ｉ＝１）を辺εとしたときに他方（ｉ＝０）を使って計算されるＬ（Ｐ）ζ（Ｐ）のメッシュ全体での総和であり、第３項は逆に辺対Ｐのうち他方（ｉ＝０）辺εとしたときにもう一方（ｉ＝１）を使って計算されるＬ（Ｐ）ζ（Ｐ）のメッシュ全体での総和である。これは、重み関数ｇ（ε）が、辺εが自身の変位に与える影響を示す重み成分と、一対の辺がそれぞれ相手の変位に影響を与えることを示す重み成分と、から構成されていることを示す。
【００９１】
【数５２】
2009301149.tif 000050
従来のように辺対Ｐを使わないで計算する場合は、四面体要素τ０で辺εの辺ベクトル
2009301149.tif 000051
（外字２）と力ｆを確定し、エレメント・バイ・エレメントで順次隣接する四面体要素τ１の計算を進めて各辺εの辺ベクトルｐ〜（ε）と力ｆを確定し、これをメッシュＭ全体で四面体要素ごとに続けて計算することになる。
【００９２】
しかし、本発明のように辺対Ｐという新たなプリミティブを採用すると、各プリミティブ、とくに辺対Ｐごとに重み関数ｇ（ε）を計算することができ、高速に演算することができる。
【００９３】
以上説明した方法を計算のフローで説明する。図３は本発明の実施の形態１の非線形の有限要素法による構造解析方法で力ｆを計算するときのフローチャートである。
【００９４】
最初に、複数の四面体要素τで分割された対象物モデルの節点νの位置ベクトルｐ（ν）を設定し、これによりｐ〜（ε）＝ｐ（νｅ（ε，０））−ｐ（νｅ（ε，１））を計算して、各四面体要素τの辺εの辺ベクトル
2009301149.tif 000052
（外字２）を決定する（ｓｔｅｐ１）。なお、ｓｔｅｐ１はｓｔｅｐ２〜７と分離された処理であり、ｓｔｅｐ２〜７と並行して行え、この前後を含んでどこで計算するのでもよい。すなわち、ｓｔｅｐ７の計算の後に処理するのでもよい。
【００９５】
次に、各四面体要素τにおいて辺εａ，辺εｂの組み合わせを指定する辺対Ｐを設定する（ｓｔｅｐ２）。その際、この四面体要素τの前に他の四面体要素τ’においてその組み合わせが既に設定されているものは除きながら、辺対Ｐを設定する。
【００９６】
予め四面体要素τの任意の辺εａ，辺εｂをパラメータとして歪みエネルギーＷτを因数分解した形で計算しておき、この因数となる関数Ｌ（εａ，εｂ）を参照し、Ｌ（ε）＝Ｌ（ε,ε）とおいて関数Ｌ（ε）を計算し、また、Ｌ（Ｐ）＝Ｌ（Ｐ,Ｐ）として関数Ｌ（Ｐ）を計算する（ｓｔｅｐ３）。
【００９７】
続いて、この歪みエネルギーＷτを構成する因数の関数ζ（ε）を参照し、関数Ｌ（ε）とこの関数ζ（ε）とによってＬ（ε）ζ（ε）を計算する（ｓｔｅｐ４）。
【００９８】
さらに、辺対Ｐで指定された辺εP（Ｐ，０），εP（Ｐ，１）ごとに関数ζ（εP（Ｐ，０）），ζ（εP（Ｐ，１））をそれぞれ計算し（ｓｔｅｐ５）、さらに、関数Ｌ（Ｐ）と関数ζ（εP（Ｐ，０）），ζ（εP（Ｐ，１））のそれぞれの積を計算して（ｓｔｅｐ６）、すべての辺対Ｐに関してこれらの和をとってｓｔｅｐ４のＬ（ε）ζ（ε）に加算することにより、辺対Ｐに関して一括して（まとまって）重み関数ｇ（ε）を計算する（ｓｔｅｐ７）。
【００９９】
この重み関数ｇ（ε）と辺ベクトル
2009301149.tif 000053
（外字２）の積ｐ〜（ε）ｇ（ε）を計算し、各辺εに関して、辺ベクトルｐ〜（ε）と重み関数ｇ（ε）の積を計算すると共に、各辺εの両端で共通する属性をもつ節点νについてこの属性に従って正負を付して積ｐ〜（ε）ｇ（ε）の総和をとる（ｓｔｅｐ８）。すなわち、辺ε自身によるｐ〜（ε）ｇ（ε）＝ｆｔｍｐと、同じ方向に位置するなど属性が共通する節点νｅ（ε，０），νｅ（ε，１）ごとのｐ〜（ε）ｇ（ε）＝ｆ（νｅ（ε，０）），ｐ〜（ε）ｇ（ε）＝ｆ（νｅ（ε，１））を属性に従って総和を計算することにより、節点νに作用する力ｆ（ν）を求めることができる。
【０１００】
これにより、グリーン・ラグランジュ歪みに基づいた非線形の有限要素法によって構造解析するときに、エレメント・バイ・エレメントの計算において発生する重複する演算を辺対Ｐという概念で無くすことができ、対象物モデルの各節点に働く力を計算するときの計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【０１０１】
図４は本発明の実施の形態１における力を計算するアルゴリズム｛ａｌｇＦ［ｐ］｝である。（数５０）と（数５２）よれば、ｆ（ν）（ν∈Ｍ）の要素から構成される力ｆは、ｐ（ν）（ν∈Ｍ）の要素から構成される位置ベクトルｐが与えられたとき、図４のアルゴリズムで計算できる。
【０１０２】
まず、メッシュMに含まれる辺εに対して、（数３５）、（数４２）、（数４８）によって１つの辺εの一端である第１の節点νｅ（ε，０）、他端である第２の節点νｅ（ε，１）間の辺ベクトル
2009301149.tif 000054
（外字２）を計算し、（数４２）によって辺εのζ（ε）を求める。これにより、（数２６）によって辺ε自身の変形による影響を示す重み関数ｇ（ε）の第１項を計算する。
【０１０３】
第１項を計算するサブルーチンを終了したら、次に（数５２）の第２項、第３項を計算するサブルーチンに進む。これらは非対角成分となる辺対Ｐに依存するものとなる。メッシュMに含まれる辺対Ｐの辺εＰのうち、一方εＰ（Ｐ，１）（すなわちｉ＝１の辺）について他方εＰ（Ｐ，０）（すなわちｉ＝０の辺）を使ってζ（εＰ（Ｐ，０））を計算し、これを変数ｇ（εＰ（Ｐ，１））に代入して初期値（重み関数ｇ（ε）の第１項）に加え、メッシュＭ全体の辺対Ｐについてこのζ（εＰ（Ｐ，０））の計算を繰り返し、順次和をとって第２項を計算する。
【０１０４】
同様に、この辺εＰの他方のεＰ（Ｐ，０）についてももう一方εＰ（Ｐ，１）を使ってζ（εＰ（Ｐ，１））を計算し、これを変数ｇ（εＰ（Ｐ，０））に代入して初期値に加え、メッシュＭ全体で辺対Ｐについてこのζ（εＰ（Ｐ，１））の計算を繰り返し、順次和をとって第３項を計算し、辺ベクトル
2009301149.tif 000055
（外字２）に対する重み関数ｇ（ε）を計算する。
【０１０５】
このサブルーチンが終了したら、以上の計算結果
2009301149.tif 000056
（外字２）、ｇ（ε）を使って力ｆ（ν）を計算する。変数ｆｔｍｐをこれらの積として（数５３）のように定義し、これを初期値としてある辺εの第１の節点νｅ（ε，０）を使ってｆ（νｅ（ε，０））を計算し、初期値に加算してこれをあらためて変数ｆｔｍｐとし、節点νｅについてｆ（νｅ（ε，０））の計算を繰り返し、ｆ（νｅ（ε，０））＋＝ｆｔｍｐの計算を行う。第２の節点νｅ（ε，１）についてもｆ（νｅ（ε，１））を計算し、初期値ｆｔｍｐに加算してこれをあらためて変数ｆｔｍｐとし、節点νｅについてｆ（νｅ（ε，１））の計算を繰り返し、ｆ（νｅ（ε，１））−＝ｆｔｍｐの計算を行う。これらから各節点νにおけるｆ（ν）を要素とする力ｆが得られる。ｆと、このときの辺ベクトル
2009301149.tif 000057
（外字２）、ｇ（ε）の値が求まったら計算を終了する。
【０１０６】
【数５３】
2009301149.tif 000058
さて、以上説明したプログラムのアルゴリズム｛ａｌｇＦ｝は、辺εの数をＮε、辺対Ｐの数をＮＰとすると、掛算の浮動点の数ＮＭＬＴ，ｆがＮＭＬＴ，ｆ＝７Ｎε＋２ＮＰ≒２６．２７Ｎτであり、加算の浮動点の数ＮＡＤＤ，ｆがＮＡＤＤ，ｆ＝１２Ｎε＋２ＮＰ≒３２．１８Ｎτとなる。ここで、プリミティブに関する節点、辺、面、四面体要素の数をそれぞれＮν、Ｎε、ＮＳ、ＮＰ、Ｎτとすると、Ｎν：Ｎε：ＮＳ：ＮＰ：Ｎτは１：６．５：１１：５．５となる。ＮＰ＝３ＮＳ＋３Ｎτであるから、Ｎν：Ｎε：Ｎτ：ＮＰ≒０．１８：１．１８：１：９程度となる。
【０１０７】
ここで、比較のため、定数行列の剛性行列Ｋを備えたｆ＝Ｋ（ｐ−ｐｉｎｉ）で表される線形の有限要素法を考えると、掛算の浮動点の数ＮＭＬＴ，ｆＬはＮＭＬＴ，ｆＬ＝９Ｎν＋１８ＮＰ≒２２．９１Ｎτであり、加算の浮動点の数ＮＡＤＤ，ｆＬはＮＡＤＤ，ｆＬ＝６Ｎε＋１８ＮＰ≒２２．３６Ｎτになる。従って、実施の形態１の計算は非線形の計算ではあるが、線形と比較して４５％程度の計算量の増加に止まっているのが分かる。
【０１０８】
これに対し、非特許文献２において報告された従来の非線形の有限要素法の計算では、２次元の計算であっても、線形の有限要素法と比較して５倍の計算量が必要である。本発明のプログラムがきわめて高速であることが分かる。
【０１０９】
以上のことから、本発明の実施の形態１のプログラムは次の２つの特徴を有していることが分かる。１つ目は、非線形のグリーン・ラグランジュ歪みを前提にしているにも関わらず、四面体要素τの歪みエネルギーＷτを（数４４）の積形式で記述して演算し易い表現にしていることである。２つ目は、エレメント・バイ・エレメントで順次隣接する四面体要素の計算を進めていく方法をやめて、辺対Ｐ（新たなプリミティブ）を使った方法にしたことである。何れの点も従来の処理を一変させるものである。そして、この場合、Ｎε＝６Ｎτ、ＮＰ≒１５Ｎτであるから、ＮＭＬＴ，ｆ＝７２Ｎτ、ＮＡＤＤ，ｆ＝１０２Ｎτとなり、計算上の浮動点の数を３７％以下に大幅に低減できることが分かる。
【０１１０】
［対象物モデル全体の接線剛性行列Ｋの計算］
続いて、（数２７）を基に対象物モデル全体の接線剛性行列Ｋについて計算する。Ｋは３Ｎν×３Ｎνの対称行列である。この３×３のブロックは次の（数５４）のように表される。ここで、ＫＰ（Ｐ）は（数５５）で定義される。
【０１１１】
【数５４】
2009301149.tif 000059
【０１１２】
【数５５】
2009301149.tif 000060
（数５４）の偏微分はνａとνｂが共通の辺εに属しているとき「０」でないから、（数５６）が得られる。
【０１１３】
【数５６】
2009301149.tif 000061
Ｋ（ν）はν∈ＭであってＫの対角行列であり、Ｋ（ε）はε∈ＭであってＫの非対角行列である。Ｋは両部分行列の和として表現されることになる。Ｋ（ν）とＫ（ε）には（数５７）の関係があるから、（数５８）を導くことができる。この中でｔ（ν，ε）は（数５９）である。
【０１１４】
【数５７】
2009301149.tif 000062
【０１１５】
【数５８】
2009301149.tif 000063
【０１１６】
【数５９】
2009301149.tif 000064
（数５８）の第３項は（数６０）と書き換えることができ、（数５８）の第４項は（数６１）と書き換えることができる。
【０１１７】
【数６０】
2009301149.tif 000065
【０１１８】
【数６１】
2009301149.tif 000066
この条件ｃｏｎｄＡは（数６２）であり、ｃｏｎｄＢは（数６３）と定義される。この２つの条件を満たすとき、εは節点νｅ（εＰ（Ｐ，０），ｉ）と節点νｅ（εＰ（Ｐ，１），ｊ）を結んだ辺となる。この辺εの向きは（数６４）（数６５）のように書き直すことができる。ここにεν（νａ，νｂ）は節点νａ，νｂが存在するときのνａ，νｂを結んだ辺である。条件ｃｏｎｄＡ（ε，Ｐ，ｉ，ｊ），ｃｏｎｄＢ（ε，Ｐ，ｉ，ｊ）の真偽値を評価して、これらの条件と重み関数ｇ（ε）、辺ベクトルｐ〜（ε）、関数Ｌ（ε）を使って（数５５），（数５８）により対象物モデル全体の接線剛性行列の部分行列Ｋ（ε）を計算する。
【０１１９】
【数６２】
2009301149.tif 000067
【０１２０】
【数６３】
2009301149.tif 000068
【０１２１】
【数６４】
2009301149.tif 000069
【０１２２】
【数６５】
2009301149.tif 000070
このｓ（Ｐ，ｉ，ｊ）は辺の向きを示すフラグで、εＰ２（Ｐ，ｉ，ｊ）がεＰ（Ｐ，１）とεＰ（Ｐ，０）のうちεＰ（Ｐ，０）側に向いていたら「１」であり、逆に向いていれば「−１」、εＰ２（Ｐ，ｉ，ｊ）が存在しない場合は「０」である。これらを使って書き直すと、（数６６）（数６７）となる。
【０１２３】
【数６６】
2009301149.tif 000071
【０１２４】
【数６７】
2009301149.tif 000072
図５は（数６４）と（数６５）の定義を示すもので、辺対Ｐが共通の節点νをもつ辺であるかないかで２通りに分かれることを示している。図５の｛タイプ−１｝の辺対Ｐは、νｅ（εＰ（Ｐ，０），１）＝νｅ（εＰ（Ｐ，１），１）であるため、νｅ（εＰ（Ｐ，１），１）が存在しない場合の辺対Ｐである。この｛タイプ−１｝の辺対Ｐは存在しないνｅ（εＰ（Ｐ，１），１）に対応するｓ（Ｐ，ｉ，ｊ）は０となる。これに対し、｛タイプ−２｝の辺対ＰはεＰ２（Ｐ，ｉ，ｊ）がすべて存在する場合を示す。
【０１２５】
この実施の形態１における非線形の有限要素法による構造解析方法で接線剛性行列Ｋを計算するときの流れについて説明すると以下のようになる。図６は本発明の実施の形態１における非線形の有限要素法による構造解析方法で対象物モデル全体の接線剛性行列Ｋを計算するときのフローチャートである。
【０１２６】
図３で得られた関数Ｌ（Ｐ）と辺εＰ（Ｐ，０），εＰ（Ｐ，１）の辺ベクトルｐ〜（εＰ（Ｐ，０）），ｐ〜（εＰ（Ｐ，０））Ｔとから（数５５）に基づいて関数ＫＰ（Ｐ）を求める（ｓｔｅｐ１１）。
【０１２７】
次に、辺εＰ（Ｐ，０）、εＰ（Ｐ，１）の各節点に対して条件ｃｏｎｄＡ（ε，Ｐ，ｉ，ｊ）、ｃｏｎｄＢ（ε，Ｐ，ｉ，ｊ）を計算し（ｓｔｅｐ１２）、条件ｃｏｎｄＡ（ε，Ｐ，ｉ，ｊ）、ｃｏｎｄＢ（ε，Ｐ，ｉ，ｊ）ごとに整数ｉ，ｊに関する（ｉ＋ｊ）を計算して関数ＫＰ（Ｐ）の正負を計算する（ｓｔｅｐ１３）。
【０１２８】
これと図３で得られた重み関数ｇ（ε）、辺ベクトル
2009301149.tif 000073
（外字２），ｐ〜（ε）Ｔ、関数Ｌ（ε）を使って（数５８）によりモデル全体の接線剛性行列Ｋ（ε）を計算する（ｓｔｅｐ１４）。
【０１２９】
これにより、エレメント・バイ・エレメントの計算において発生する重複する演算を辺対Ｐという概念で無くすことができ、モデル全体の接線剛性行列を計算するときの計算量を低減することができる。動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【０１３０】
以上説明した（数５８）、（数６０）、（数６１）を使って、対象物モデル全体の接線剛性行列Ｋの非対角成分は、図７に示すアルゴリズム｛ａｌｇＫＥ｝で計算することができる。
【０１３１】
図７において、メッシュＭに属するεに対して、（数５８）の第１項と第２項を計算してＫ（ε）と定義する。次に、メッシュＭに属するＰに対して、（数５５）を計算してＫＰ（＝ＫＰ（Ｐ））と定義する。
【０１３２】
Ｋ（ε）の計算を行うため、辺εＰの一方の辺に存在する番号ｉの節点νｅと、辺εＰの他方の辺に存在する番号ｊの節点νｅとをつなぐ辺εＰ２（Ｐ，ｉ，ｊ）をεＰ２（Ｐ，ｉ，ｊ）＝εν（νｅ（εＰ（Ｐ，０），ｉ），νｅ（εＰ（Ｐ，１），ｊ）））と定義し、ｉ，ｊを０，１の間で変化させる。ｓ（Ｐ，ｉ，ｊ）＝１の場合は、このｉ，ｊについて（−１）ｉ＋ｊＫＰを求め、これにＫ（εＰ２（Ｐ，ｉ，ｊ））を順次加算してＫ（ε）の第３項とする。また、ｓ（Ｐ，ｉ，ｊ）＝−１の場合は、このｉ，ｊについて（−１）ｉ＋ｊＫＰＴを求め、これにＫ（εＰ２（Ｐ，ｉ，ｊ））を順次加算してＫ（ε）の第４項とする。ｓ（Ｐ，ｉ，ｊ）＝０の場合は０で、Ｋ（ε）に変化はない。
【０１３３】
なお、以上説明したアルゴリズム｛ａｌｇＫＥ｝は辺ベクトル
2009301149.tif 000074
（外字２）、ｇ（ε）を用いて演算するため、図４で示した辺ベクトル
2009301149.tif 000075
（外字２）、ｇ（ε）を計算する第１のアルゴリズム｛ａｌｇＦ［ｐ］｝を計算した後で計算する必要がある。
【０１３４】
次に、対象物モデル全体の接線剛性行列Ｋの対角成分Ｋ（ν）は図８の第３のアルゴリズム｛ａｌｇＫＶ｝で計算することができる。図８に示すように、Ｋ（ν）の初期値を零行列とし、（数５７）に従って、辺εの一端となる節点νｅ（ε，０）のＫ（νｅ（ε，０））を計算する。これを初期値から引き、この結果をあらためてＫ（ν）とし、辺εについてこの計算を繰り返し、Ｋ（ν）の値を更新する。同様に、辺εの他端となる節点νｅ（ε，１）のＫ（νｅ（ε，１））を計算し、これを初期値から引き、この結果をあらためてＫ（ν）とし、辺εについてこの計算を繰り返す。これにより、Ｋ（ν）が得られる。
【０１３５】
ところで、Ｋ（ν）は対称マトリクスであって、６個の独立な要素しか有していない。このとき加算の浮動点の数は１２Nｅ（辺εの数Nｅ）であり、第２，第３のアルゴリズムによる積算の浮動点の数ＮＭＬＴ，Ｋと加算の浮動点の数ＮＡＤＤ，Ｋは、それぞれＮＭＬＴ，Ｋ＝９Ｎε＋１２ＮＰ≒１１８．６４Ｎτ、ＮＡＤＤ，Ｋ＝１５Ｎε＋２７ＮＰ１＋３６ＮＰ２≒２８７．７３Ｎτと評価できる。
【０１３６】
辺εと辺対Ｐを使った場合の積算と加算と浮動点の数と、従来のエレメント・バイ・エレメントによる積算と加算と浮動点の数を比較するには、１個の四面体要素τに対して従来のエレメント・バイ・エレメントでは６個の辺と１５個の辺対があることを考慮して、Ｎεに６Ｎτ、ＮＰに１５Ｎτを代入して比較すればよい。すると、積算の浮動点の数ＮＭＬＴ，Ｋは２３４Ｎτ、加算の浮動点の数ＮＡＤＤ，Ｋは５２２Ｎτとなり、計算量が５８％に減ることが分かる。
【０１３７】
（表１）は力ｆを計算するとき、線形（ｌｉｎｅａｒ）による計算、辺対Ｐを使った本発明による計算、本発明をエレメント・バイ・エレメント（τｂτ）方式に施した計算、エレメント・バイ・エレメント（τｂτ）による計算により、１個の四面体要素当たりに必要な各計算量を示す。これをみると、エレメント・バイ・エレメント（τｂτ）による全計算量と比較して、辺対Ｐによる全計算量は１／６．６２倍であって、格段に計算量が低下している。
【０１３８】
【表１】
2009301149.tif 000076
また、（表２）は力ｆと対象物モデル全体の接線剛性行列Ｋを計算するとき、辺対Ｐを使った本発明による計算、本発明をエレメント・バイ・エレメント（τｂτ）方式に施した計算、エレメント・バイ・エレメント（τｂτ）による計算により、１個の四面体要素τ当たりで必要な演算量を示す。これをみると、エレメント・バイ・エレメント（τｂτ）による全計算量と比較して、辺対Ｐによる方法の全計算量は１／３．７６倍であって、各段に計算量が低下している。
【０１３９】
【表２】
2009301149.tif 000077
［体積歪みエネルギーの計算］
グリーン・ラグランジュ歪みのテンソルは、回転した場合だけでなく、反転するときにも不変性を保つ。このときモデルが圧縮された形状となってしまうので、これを防いで物理的な矛盾をなくすため、（数６８）のように歪みエネルギー密度ｗに付加的な歪みエネルギーΨｖｏｌ（θ）を加える。θは体積歪みであってθ＝（ｄｅｔ（Ｆ）−１）／６である。実施の形態１においては、Ψｖｏｌ（θ）としてΨｖｏｌ（θ）＝Ｋｖｏｌ｜ｍｉｎ（０，θ）｜３／３を選択している。ここにＫｖｏｌ＞０とする。さらにνをポアソン比としてＫｖｏｌ＝Ｅ／（４（１−２ν））とした。
【０１４０】
【数６８】
2009301149.tif 000078
この付加的な歪みエネルギーΨｖｏｌ（θ）が節点と接線剛性行列Ｋに与える作用について説明する。四面体要素τの付加的な歪みエネルギーＷτｖｏｌ（θ）は（数６９）で表される。ここでΨτｖｏｌ（ｃ）は（数７０）であり、Ｃτは四面体要素τの体積の６倍である。
【０１４１】
【数６９】
2009301149.tif 000079
【０１４２】
【数７０】
2009301149.tif 000080
このとき節点νに作用する力ｆτｖｏｌ（ν）は（数７１）で表され、接線剛性行列Ｋτｖｏｌ（ε）は（数７２）で表される。この中でΔＣは、ΔＣ＝Ｃτ−Ｃτｉｎｉである。
【０１４３】
【数７１】
2009301149.tif 000081
【０１４４】
【数７２】
2009301149.tif 000082
辺εと節点νの方向と添え字を図２の四面体要素τ０のように定義すると、（数７３）（数７４）のように表現でき、対象物モデル全体の力ｆτｖｏｌ（ν）、接線剛性行列Ｋτｖｏｌ（ε）は図９の第４のアルゴリズム｛ａｌｇＢτ｝で計算できる。
【０１４５】
【数７３】
2009301149.tif 000083
【０１４６】
【数７４】
2009301149.tif 000084
その計算の手順は、アルゴリズム｛ａｌｇＢτ｝の中で添え字「＊」は適切な整数とし、ε＊の
2009301149.tif 000085
（外字２）を使って図９のようにｃ０，ΔＣτ，Ψ１，Ψ２を定義する。ｃ０と同様にｃｊ（ｉ＝１，２）を計算し、（数７３）とｃｊ（ｉ＝０，１，２）からｆτｖｏｌ（νｊ）を求めると共に、ｆτｖｏｌ（ν３）＝−ｆτｖｏｌ（ν０）−ｆτｖｏｌ（ν１）−ｆτｖｏｌ（ν２）を計算する。そして、ｊ＝０，１，２，３とすると、ｆ（νｊ）はｆτｖｏｌ（ν０）にｊ＝０，・・，（ｊ−１）までのｆτｖｏｌ（νｊ）を加算したｆ＝Σｊ＝０，１，・・，（ｊ−１）ｆτｖｏｌ（νｊ）になる。これにより対象物モデル全体の力ｆが求まる。
【０１４７】
続いて、ｉ＝０，１，２と、ｉ＝３，４，５に分けてｘｉを計算し、ｙｉ（ｉ＝０，１，２，３）を計算する。これを基にｃ＊ｙ＊Ｔ＋（−１）＊［ｘ＊×］を計算し、ｉ＝０，１，・・，（ｉ−１）に対してＫ（εｉ）を順次加算し、Ｋ（εｉ）をΣｉ＝０，１，・・，（ｉ−１）Ｋ（εｉ）として算出する。［ｘ＊×］はｘ＊の外積計算を表す行列を示す。これにより対象物モデル全体の接線剛性行列Ｋ（ε）が求まる。
【０１４８】
ここで、第４のアルゴリズム｛ａｌｇＢ｝τは（１０２＋ｎｍψ）個の積算回数と、（１３５＋ｎａψ）個の加算回数を含んでいる。なお、このｎｍψは
2009301149.tif 000086
（外字５）の積算の回数であり、ｎａψは
2009301149.tif 000087
（外字６）の積算の回数である。アルゴリズム｛ａｌｇＢτ｝の計算量と表２の全計算量との比をとると、概ね（２３７＋ｎｍψ＋ｎａψ）／４６４．８２＜０．５５であって、エレメント・バイ・エレメントで計算した全計算量の３．７６倍に比較すれば、大幅に計算量を低減できる。
【０１４９】
このアルゴリズム｛ａｌｇＢτ｝は、アルゴリズム｛ａｌｇＫＥ｝の後で、しかもアルゴリズム｛ａｌｇＫＶ｝の前に実行しなければならない。従って、節点に作用する力ｆと、対象物モデル全体の接線剛性行列Ｋを演算する全アルゴリズム｛ａｌｇＡＬＬ｝は図１０に示すようになる。第１のアルゴリズム｛ａｌｇＦ［ｐ］｝を実行し、次に第２のアルゴリズム｛ａｌｇＫＥ｝を実行する。続いて、第３のアルゴリズム｛ａｌｇＢτ｝を実行してから第４の｛ａｌｇＫＶ｝を実行する。これによってｆとＫとが演算できる。
【０１５０】
次に、実施の形態１におけるシミュレーション装置について説明する。図１１はグリーン・ラグランジュ歪みに基づいた非線形の有限要素法によって構造解析を行う実施の形態１におけるシミュレーション装置である。本発明にいうシミュレーション装置とは、シミュレーションを行ういわゆるシミュレータだけでなく、３Ｄモデル教育やコンピュータゲーム、あるいは実際の動作に近いアニメーションの作成など、モデルを動作させるインタラクティブなアプリケーションを搭載する装置すべてを包含するものである。
【０１５１】
図１１において、実施の形態１のシミュレーション装置は、少なくとも上述した力ｆ及び／又は対象物モデル全体の接線剛性行列Ｋをコンピュータに計算させるプログラムが格納されたハードディスクやフラッシュメモリ等の記憶手段１００と、触覚入力デバイスなどのモデルに与える変形条件や外力などの量を入力できる入力手段２００と、ハードウェアとして中央演算処理装置（ＣＰＵ）を備え、ＲＡＭ（Read Random Memory）に上記プログラムを含むプログラムを読み出して、逐次演算を実行する機能実行手段としての制御手段３００と、入力手段２００からの入力に応じて計算結果の変形状態を表示する表示手段４００とから構成されている。なお、入力手段２００からの入力は、図示しないインターフェースにおいて感覚としての入力と演算値とのマッチングが図られる。
【０１５２】
このグリーン・ラグランジュ歪みに基づいて構造解析を行うシミュレーション装置では、エレメント・バイ・エレメントの計算において発生する重複する演算を辺対Ｐという概念により無くすことができ、非線形の有限要素法による計算量を低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。この構成によって、入力手段２００からモデルに変位を加えたとき、表示手段４００にモデルが変形される状態を実時間で表示することができる。
【０１５３】
（実施例１）
具体的にＮν個の節点を有するメッシュＭで運動方程式（数７５）を解いた結果を説明する。所定の質量をもった物体が外力と抵抗を受けながら変形する一般式である。
【０１５４】
【数７５】
2009301149.tif 000088
ｆは変形に伴う物体内部の力ｆ＝ｆ（ｐ）∈Ｒ３Ｎνであり、Ｍは質量、Ｂは粘性係数、ｈ∈Ｒ３Ｎνは節点に作用する外力の成分である。レイリー減衰モデルを採用してＢをＢ＝αＭ＋βＫとした。α，βは定数であり、質量Ｍは対角行列とし、ＫはＫ＝Ｋ（ｐ）であて対象物モデル全体の接線剛性行列である。外力ｈは変形条件あるいは触覚入力デバイスとの間をスプリングダンパーで接続したときに作用する力と仮定し、外力ｈはデバイスの状態（位置または速度）ｑ、及び位置ｐによって決定されｈ＝Ｈ（ｐ，ｖ，ｑ）となる。
【０１５５】
ここでＨｐ＝∂Ｈ／∂ｐ、Ｈｖ＝∂Ｈ／∂ｖとおき、（数７５）の微分方程式を図１２に従って後退（陰的）オイラー法で数値計算する。図１２の中でＡｖ（ｋ）＝ｂはＰＣＧＭ法（preconditioned conjugate gradient method）で解いている。なお、Ａはほとんどの場合において正定行列であり、正定行列であるという仮定の下でＰＣＧＭ法を採用した。図１２のＴはこのときのステップ時間の長さである。
【０１５６】
（表３）に記載したように板状モデル、狼状モデル、兎状モデルの３種類のモデルのメッシュを使って計算を行った。メッシュに対するヤング率とポアソン比は０．０５ＭＰａ、０．４とした。Ｍは対角要素が密度ρ＝１０−６Ｋｇ／ｍｍ３が分布しているものと近似した。αはα＝０．５ｓ−１、βはβ＝０．１、ＴはＴ＝３０ｍｓとして好適な値を選択した。
【０１５７】
【表３】
2009301149.tif 000089
図１３（ａ）は板状モデルの初期の状態を示しており、図１３（ｂ）はこれを大きく変形して折り曲げた状態を示している。また、図１３（ｃ）は板状モデルに大きなねじりを加えた状態を示している。折り曲げた状態、大きなねじりを加えた状態の何れにおいても、従来の線形の有限要素法のモデルでは、歪みが大きくなると体積が増加して不自然な形状になっていたが、実施例１で計算した板状物体は大きな変位、大きな回転であるにも関わらず、こうした異常がみられず、物理的にまた実感として自然な変形となっている。しかも、非線形の有限要素法の計算としては高速になる。狼状物体、兎状物体についても、図示はしないが、大きな変位、大きな回転に対して異常がみられず、物理的にまた実感として自然な変形となった。
【０１５８】
このように本発明の非線形の有限要素法による構造解析方法、プログラム、記録媒体、シミュレーション装置によれば、非線形の有限要素法による計算量を大幅に低減することができ、動的な計算においても３次元モデルの実時間の動作を可能にすることができる。
【産業上の利用可能性】
【０１５９】
本発明は、非線形の有限要素法による構造解析のプログラム、シミュレーション装置に適用できる。
【図面の簡単な説明】
【０１６０】
【図１】本発明の実施の形態１における非線形の有限要素法の構造解析方法が前提とする対象物モデルのメッシュの説明図
【図２】四面体要素の一面を共通にして接続した各プリミティブの説明図
【図３】本発明の実施の形態１の非線形の有限要素法による構造解析方法で力ｆを計算するときのフローチャート
【図４】本発明の実施の形態１におけるプログラムの各節点の力を演算するアルゴリズムの説明図
【図５】辺対が共通の節点をもつ辺か否かで２通りのタイプに分かれることの説明図
【図６】本発明の実施の形態１における非線形の有限要素法による構造解析方法で対象物モデル全体の接線剛性行列を計算するときのフローチャート
【図７】本発明の実施の形態１におけるプログラムの接線剛性行列の非対角成分を演算するアルゴリズムの説明図
【図８】本発明の実施の形態１におけるプログラムの接線剛性行列の対角成分を演算するアルゴリズムの説明図
【図９】本発明の実施の形態１におけるプログラムの節点に作用する力と対象物モデル全体の接線剛性行列を演算するアルゴリズムの説明図
【図１０】本発明の実施の形態１におけるプログラムの一連のアルゴリズムを演算するアルゴリズムの説明図
【図１１】本発明の実施の形態１におけるシミュレーション装置の構成図
【図１２】本発明の実施の形態１におけるプログラムの後退（陰的）オイラー法で数値計算を行うアルゴリズムの説明図
【図１３】（ａ）板状モデルの初期の状態を示す図、（ｂ）板状モデルを大きく変形して折り曲げた状態を示す図、（ｃ）板状モデルに大きなねじりを加えた状態を示す図
【符号の説明】
【０１６１】
Ｍ    メッシュ
τ    四面体要素
Ｓ    面
ε，εａ，εｂ    辺
ν    節点
ｆ    力
ｐ（ν）    位置ベクトル
Ｋ（ｐ）    接線剛性行列
１００    記憶手段
２００    入力手段
３００    制御手段
４００    表示手段
【図１】
2009301149.tif 000104
【図２】
2009301149.tif 000105
【図３】
2009301149.tif 000106
【図４】
2009301149.tif 000107
【図５】
2009301149.tif 000108
【図６】
2009301149.tif 000109
【図７】
2009301149.tif 000110
【図８】
2009301149.tif 000111
【図９】
2009301149.tif 000112
【図１０】
2009301149.tif 000113
【図１１】
2009301149.tif 000114
【図１２】
2009301149.tif 000115
【図１３】
2009301149.tif 000116