(57)【要約】【課題】局所座標系に関連する、システムの一節点における少なくとも１つの解放された自由度を伴う、予め定められた荷重および変位境界条件を受けたシステムの挙動を判定する。【解決手段】プロセスは有限要素解析を用いる。幾何学的非線形解法スキームでは、プロセスは複数のインクリメンタルステップを有する。各インクリメンタルステップでは、システムの有限要素に関し、有限要素の節点の変位ベクトルを計算する。この計算は、予め定められた荷重と予め定められた変位の一部とを受けたときのシステムのエネルギを最小化することにより実行される。幾何学的線形解法スキームでは、プロセスは１回のソルビングステップにより実行される。このステップでは、システムの有限要素に関し、有限要素の節点の変位ベクトルを計算する。この計算は、予め定められた荷重と予め定められた変位とを受けたときのシステムのエネルギを最小化することにより実行される。【選択図】図６
2006294038.tif 000002
【特許請求の範囲】
【請求項１】
有限要素解析を用いて、予め定められた荷重および変位境界条件を受けたシステムの挙動を判定するためのプロセスであって、前記予め定められた変位境界条件は、前記システムの一節点における少なくとも１つの自由度の解放を含み、前記プロセスは、複数のインクリメンタルステップを備え、各インクリメンタルステップは、予め定められた変位の一部を受けたときの前記システムのエネルギを最小化することによって、前記システムの有限要素に関して、前記有限要素の節点における変位ベクトルを計算することを備え、前記最小化するステップは、解放された自由度を有する節点に関して、ローカルリファレンスフレームを用いて実行され、前記節点における解放された自由度は、前記ローカルリファレンスフレームの１つの軸に沿って表現されることを特徴とするプロセス。
【請求項２】
有限要素解析を用いて、予め定められた荷重および変位境界条件を受けたシステムの挙動を判定するためのプロセスであって、前記予め定められた変位境界条件は、前記システムの一節点における少なくとも１つの自由度の解放を含み、前記プロセスは、予め定められた変位を受けたときの前記システムのエネルギを最小化することによって、前記システムの有限要素に関して、前記有限要素の節点における変位ベクトルを計算するステップを備え、前記最小化するステップは、解放された自由度を有する節点に関して、ローカルリファレンスフレームを用いて実行され、前記節点における解放された自由度は、前記ローカルリファレンスフレームの１つの軸に沿って表現されることを特徴とするプロセス。
【請求項３】
前記システムは、前記変位ベクトルが並進コンポーネントのみを有する連続要素を用いてモデル化されることを特徴とする請求項１または２に記載のプロセス。
【請求項４】
前記システムは、前記変位ベクトルが並進コンポーネントおよび回転コンポーネントを有する構造要素を用いてモデル化されることを特徴とする請求項１または２に記載のプロセス。
【請求項５】
節点リファレンスフレーム（４３，４８）が有限要素（４０）の節点（４１，４６）に提供され、前記計算するステップは、前記節点（４１，４６）における前記変位ベクトルの回転コンポーネント分だけ前記節点リファレンスフレーム（４３，４８）を回転することによって、前記有限要素の変位した節点（７１，７６）に関して、最新の節点リファレンスフレーム（７３，７８）を計算することをさらに備えることを特徴とする請求項１ないし４のいずれかに記載のプロセス。
【請求項６】
前記有限要素に関する前記節点の位置と、前記節点リファレンスフレームとに基づいて、前記システムに合わせたスムーズ形状を補間するステップをさらに備えることを特徴とする請求項５に記載のプロセス。
【請求項７】
補間された前記スムーズ形状を表示するステップをさらに備えることを特徴とする請求項６に記載のプロセス。
【請求項８】
前記システムの前記最小化エネルギは、節点における節点リファレンスフレームと、前記節点を含む有限要素のエレメンタリリファレンスフレームとの差の関数であることを特徴とする請求項１ないし７のいずれかに記載のプロセス。
【請求項９】
前記計算するステップの間、エレメンタリリファレンスフレームを節点リファレンスフレームに変換する、有限要素の節点における回転を計算することをさらに備え、前記最小化エネルギは、計算された前記回転に依拠することを特徴とする請求項８に記載のプロセス。
【請求項１０】
前記システムの前記最小化エネルギは、節点における節点リファレンスフレームと前記節点を含む有限要素のエレメンタリリファレンスフレームとの差と、スケーリング係数との積の関数であり、前記プロセスは、一インクリメンタルステップと、別のインクリメンタルステップとで前記スケーリング係数を変更することをさらに備えることを特徴とする請求項８または９に記載のプロセス。
【請求項１１】
前記スケーリング係数は手動で変更されることを特徴とする請求項１０に記載のプロセス。
【請求項１２】
前記スケーリング係数は、インクリメンタルステップにおいて前記エネルギを最小化するために、反復回数に応じて変更されることを特徴とする請求項１０に記載のプロセス。
【請求項１３】
前記スケーリング係数は、有限要素用のエレメンタリリファレンスフレームと、前記有限要素に関する節点リファレンスフレームとの差に応じて変更されることを特徴とする請求項１０に記載のプロセス。
【請求項１４】
前記エネルギを前記最小化するステップは反復的に実行され、該ステップは、
余剰変動と、変位変動とのスカラ積を計算することと、
前記スカラ積の値に応じて、反復を中止することとを備えることを特徴とする請求項１ないし１３のいずれかに記載のプロセス。
【請求項１５】
請求項１ないし１４のいずれかに記載の、コンピュータで実現されるプロセスの全ステップを実行するのに適したコンピュータプログラムコード手段を備えた、コンピュータ読取り可能媒体上に存在するコンピュータプログラムの製品。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピュータプログラムおよびコンピュータシステムの分野に関し、より詳細には、部品設計プログラム、および、そのような制約ベースのシステム（constraints-based system）に関する。
【背景技術】
【０００２】
部品または部品のアセンブリの設計用に、本件出願人がＣＡＴＩＡ（登録商標）という商標の下で提供しているものなど、多数のシステムおよびプログラムが市販されている。これらのいわゆるコンピュータ支援設計（ＣＡＤ）システムによって、ユーザは、部品または部品のアセンブリの複雑な３次元（３Ｄ）モデルを構築し、操作することができる。これらのシステムおよびプログラムは、モデルを定義するために様々な制約を用いる。モデルが編集されるとき、制約のセットはシステムによって解かれる。制約のセットを解くために使用されるプログラムまたはシステムは、一般に「ソルバ」と呼ばれる。このようなソルバは、ＣＡＤ／ＣＡＭ／ＣＡＥシステム、または、より一般的には、任意の種類のオブジェクト（object）を定義するために制約を用いる任意のシステムにおいて、使用される。ソルバは、設計されているオブジェクトと、こうしたオブジェクトに与える種類の制約とに適合する。ＣＡＴＩＡ（登録商標）において使用されるものなどのソルバは、寸法制約、圧迫、オブジェクト間の接触などを含む制約を伴うソリッドオブジェクトの設計用に適している。
【０００３】
予め定められた荷重および変位境界条件を受けた、変形しやすいソリッドフィジカルシステム（deformable solid physical system）の変形した形状をシミュレートすることを可能にするソルバが必要とされる。
【０００４】
ＣＡＤ／ＣＡＭ／ＣＡＥシステムにおいて通常設計される一部のフレキシブルシステム（flexible system）、および、オブジェクトは、個々の点が荷重の作用下で任意に大きな回転を個別に受ける場合があるという点で、それらの高変形性（「幾何学的非線形性」）によって特徴付けられる。すなわち、そのようなフレキシブルシステムでは、システムの変形した形状と、システムの変形していない形状とが大幅に異なり得る。それに対し、既存のＣＡＤ／ＣＡＭ／ＣＡＥシステムにおける幾何学的線形ソルバは、オブジェクトの変形した形状と、変形していない形状とがほとんど変わらないという仮定の下で、正常に機能する。
【０００５】
より一般的には、タンジェントによりラジアンを単位とする回転角度の近似値を求めることができないような大回転（通常は５度であるが、この数字は、許容されるモデリングエラー量に応じて変わる可能性がある）を受けるシステムは、通常「幾何学的非線形システム」とみなされる。あるいは、グリーン−ラグランジュ測定（Green-Lagrange measure）などの「大ひずみ」定義
【０００６】
【数１】
2006294038.tif 000003
【０００７】
を、ひずみテンソルコンポーネントの「微小ひずみ」定義
【０００８】
【数２】
2006294038.tif 000004
【０００９】
で置き換えることができないとき、すなわち、グリーン−ラグランジュ測定における積項
【００１０】
【数３】
2006294038.tif 000005
【００１１】
を無視することができないときに、システムは非線形であると言うことができる。
【００１２】
一般に、ソルバにおいて、システムのいくつかの制約された（restrained）点における、自由度の部分的な解放（partial release of degrees of freedom）を可能にする必要性が存在する。例えば、システムには、すべり点（sliding point）が提供され得る。そのすべり点において、システムの移動は２方向に制約されるが、システムは、１つの方向に沿ってすべることができ、おそらくはすべり方向のまわりを回転することができる。この例では、システムは、そのすべり点において、並進について１自由度を有し、回転について１自由度を有する。以下で詳細に説明するように、従来技術（例えば、「Use of the FEM for the Design of Flexible Parts, C.A. de Hillerin, Proceedings of NAFEMS World Congress 1999, pp. 345-356」）では、そのような状況を考慮するために、制約をファクトライズ（factorise）することが教示されている。
【００１３】
非線形システムであるスレンダボディフレキシブルシステム（slender-body flexible system）に関して存在するこの必要性は、システムの変形した形状と、変形していない形状とがほとんど変わらないとみなされる線形システムについても生じる。
【００１４】
スレンダボディフレキシブルシステムの別の特徴は、有限要素解析におけるシステムの要素が、並進コンポーネントおよび回転コンポーネントを有する変位を受けるということである。この要素は、いわゆる「構造」要素（structural element）である。構造要素の他の例には、梁（beam）およびシェル（shell）が含まれる。制約された点における自由度の部分的な解放を可能にする必要性は、システムの要素に関する変位ベクトルが並進コンポーネントのみを有する、いわゆる「連続」要素（continuous element）を用いてモデル化されたシステムにおいても存在する。連続システム（continuous system）の例には、バー（bar）、膜（membrane）およびソリッド（solid）が含まれる。構造システム（structural system）と連続システムとの違いは、「Finite Elements Procedures in Engineering Analysis, K-J. Bathe, Prentice-Hall 1982」において説明されている。言うまでもなく、複雑なシステムは、同一モデル内に両方のタイプの要素を含み得る。
【発明の開示】
【発明が解決しようとする課題】
【００１５】
特に、解放された自由度（released degrees of freedom）がグローバルリファレンスフレームの軸に対応しないとき、制約をファクトライズするという従来技術による解法は、複雑で時間のかかる計算を伴う場合がある。したがって、システムのいくつかの制約された点における自由度の部分的な解放を可能にすると同時に、計算をより単純にするソルバが必要とされる。
【課題を解決するための手段】
【００１６】
本発明の一実施形態にしたがうと、上記の問題に対する解決策が提供される。この解決策は、構造要素定式化（structural element formulation）および連続体要素定式化（continuum element formulation）の両方に適用することができ、非線形解法スキーム（non-linear solution scheme）および線形解法スキーム（linear solution scheme）の両方に適用することができる。
【００１７】
連続体システム（continuum system）に関して、本発明の第１の実施形態は、静的な有限要素法のフレームワーク内で、線形の制約方程式系を処理できる方法に関するものである。この場合、部分的に制約された節点のグローバルな並進自由度に関連する制約関係（constraint relation）の数は、１または２である。これらの制約関係は、グローバルな座標軸には一致しない１つまたは２つの斜方向に沿って、特定の節点の並進変位が部分的に制約されることを表す。
【００１８】
そのような状況に対して、この第１の実施形態は、標準の「制約ファクトリゼーション（constraints factorization）」プロシージャを、新たなプロシージャによって置き換えることを提案する。
【００１９】
標準の「制約ファクトリゼーション」プロシージャは、グローバルな変位コンポーネントに関して記述された一連の連結制約方程式（coupled constraint equation）を生成することと、次いで、これを、すべての依存自由度（dependent degrees of freedom）が取り除かれた等価な形式に変換することと、次いで、（アセンブリに先立って、グローバルリファレンスフレームを用いて回転されたエレメンタリ剛性行列（elementary stiffness matrix）および荷重ベクトル（load vector）から構成されるものとみなされる）グローバルな平衡方程式系において、この形式を代用することと、最後に、この系を、その一連の連結制約方程式の零空間（nullspace）に射影する（project）ことによって、すべての依存自由度と未知のラグランジュ乗数とを取り除き、グローバルな平衡方程式系を、コンパクトかつ可逆的な方程式系に変形することと、からなる。
【００２０】
第１の実施形態にしたがう新たなプロシージャでは、まず、部分的に制約された各節点のグローバルな並進変位のコンポーネントが、ローカルな斜めの制約方向（restraint direction）から構成された正規直交型のリファレンスフレーム（orthonormal reference frame）に関連するローカルコンポーネントに変換される。ローカルコンポーネントを用いると、そのような部分的に制約された節点の並進変位に関連する制約関係が、次いで、より単純な分離した（uncoupled）形式で表現される。対応する平衡方程式系を構成するには、アセンブリに先立って、「ハイブリッド（hybrid）」エレメンタリ剛性行列および荷重ベクトルの計算が必要となる。
【００２１】
「ハイブリッド」エレメンタリ剛性行列および荷重ベクトルは、「ハイブリッド」行列回転およびベクトル回転によって得られる。すなわち、「ハイブリッド」エレメンタリ剛性行列は、様々なリファレンスフレームにおける様々な節点に対応する剛性行列の対角ブロックを選択的に回転して、その結果を非対角ブロックに適切に反映することによって得られ、荷重ベクトルも類似の処理によって得られる。
【００２２】
適切な平衡方程式系が組み立てられると、組み立てられた剛性行列および荷重ベクトルの直接的な検証（direct inspection）に基づく、単純な行および列の除去操作によって、分離した制約方程式が考慮されるようになる。結果として生じる方程式系は、「ハイブリッド」と呼ばれる。というのは、その未知数は１つだけのリファレンスフレームに関連する変位コンポーネントではないからである。それでも、すべての依存自由度と未知のラグランジュ乗数とが、取り除かれており、変形した平衡方程式系は、この場合もまたコンパクトかつ可逆的である。
【００２３】
このプロセスは、非線形の解法のみならず、インクリメンタルステップがないことが、制約が１回しか考慮される必要がないことを意味する、線形の解法に対しても適用される。
【００２４】
（以下の例で説明する「ケーブル」要素に加えて、梁要素、プレート要素およびシェル要素を含む）構造要素についてみると、連続体要素との最大の違いは、構造要素が、節点の並進自由度に加えて、節点の回転自由度も表すことにある。構造要素に関して、本発明の第２の実施形態は、静的な有限要素法のフレームワーク内で、線形の制約方程式系を処理できる方法に関するものである。この場合、部分的に制約された節点のグローバルな並進自由度に関連する制約関係の数は、１または２であり、かつ／または、部分的に制約された節点のグローバルな回転自由度に関連する制約関係の数は、１または２である。これらの制約関係は、グローバルな座標軸には一致しない斜方向に沿って、特定の節点の並進変位、および／または、回転変位が部分的に制約されることを表す。
【００２５】
そのような状況に対して、本発明の第２の実施形態は、標準の「制約ファクトリゼーション」プロシージャを、上述したものに類似した新たなプロシージャによって置き換えることを提案する。特定の差異（これは、線形化した回転キネマティックス（linearized rotation kinematics）を伴う構造システムが、非線形のインクリメンタルな解法の場合には、アップデートを必要とする、という事実に関連する。）を除けば、部分的な緩和（partial relaxation）を伴う斜めの制約（oblique constraint）を処理するためのハイブリッド化プロセス（hybridation process）に関して、他のすべてが同一に保たれる。
【００２６】
具体的には、本発明は、第１の実施形態において、有限要素解析を用いて、予め定められた荷重および変位境界条件を受けたシステムの挙動を判定するためのプロセスを提供する。予め定められた変位境界条件は、システムの一節点における少なくとも１つの自由度の解放を含む。このプロセスは、複数のインクリメンタルステップを備える。各インクリメンタルステップは、予め定められた変位の一部を受けたときのシステムのエネルギを最小化することによって、システムの有限要素に関して、その有限要素の節点における変位ベクトルを計算することを備える。最小化するステップは、解放された自由度を有する節点に関して、ローカルリファレンスフレームを用いて実行され、前記節点における解放された自由度は、ローカルリファレンスフレームの１つの軸に沿って表現される。
【００２７】
本発明は、第２の実施形態において、有限要素解析を用いて、予め定められた荷重および変位境界条件を受けたシステムの挙動を判定するためのプロセスを提供する。予め定められた変位境界条件は、システムの一節点における少なくとも１つの自由度の解放を含む。このプロセスは、予め定められた変位を受けたときのシステムのエネルギを最小化することによって、システムの有限要素に関して、その有限要素の節点における変位ベクトルを計算するステップを備える。最小化するステップは、解放された自由度を有する節点に関して、ローカルリファレンスフレームを用いて実行され、前記節点における解放された自由度は、ローカルリファレンスフレームの１つの軸に沿って表現される。
【００２８】
第１の実施形態および第２の実施形態において、以下のさらなる特徴のうちの１つまたは複数が提供され得る。
【００２９】
−システムは、変位ベクトルが並進コンポーネントのみを有する連続要素を用いてモデル化される。
【００３０】
−システムは、変位ベクトルが並進コンポーネントおよび回転コンポーネントを有する構造要素を用いてモデル化される。
【００３１】
−節点リファレンスフレーム（nodal frame of reference）が有限要素の節点に提供され、計算するステップは、節点における変位ベクトルの回転コンポーネント分だけ節点リファレンスフレームを回転することによって、有限要素の変位した節点に関して、最新の節点リファレンスフレームを計算することをさらに備える。
【００３２】
−プロセスは、有限要素に関する節点の位置と、節点リファレンスフレームとに基づいて、システムに合わせたスムーズ形状（smoothed shape）を補間するステップをさらに備える。
【００３３】
−プロセスは、補間されたスムーズ形状を表示するステップをさらに備える。
【００３４】
−システムの最小化エネルギは、節点における節点リファレンスフレームと、その節点を含む有限要素のエレメンタリリファレンスフレームとの差の関数である。
【００３５】
−プロセスは、計算するステップの間、エレメンタリリファレンスフレームを節点リファレンスフレームに変換する、有限要素の節点における回転を計算することをさらに備え、最小化エネルギは、計算された回転に依拠する。
【００３６】
−システムの最小化エネルギは、節点における節点リファレンスフレームとその節点を含む有限要素のエレメンタリリファレンスフレームとの差と、スケーリング係数（scaling factor）との積の関数であり、プロセスは、一インクリメンタルステップと、別のインクリメンタルステップとでスケーリング係数を変更することをさらに備える。スケーリング係数は、手動で変更されてもよいし、インクリメンタルステップにおいてエネルギを最小化するために、反復回数に応じて変更されてもよい。または、スケーリング係数は、有限要素用のエレメンタリリファレンスフレームと、その有限要素に関する節点リファレンスフレームとの差に応じて変更されてもよい。
【００３７】
−エネルギを最小化するステップは、反復的に実行され、このステップは、余剰変動と、変位変動とのスカラ積を計算すること、および、そのスカラ積の値に応じて、反復を中止することを備える。
【００３８】
最後に、本発明は、コンピュータ上でこのようなプロセスの全ステップを実行するのに適したコンピュータプログラムコード手段を備える、コンピュータ読取り可能媒体上に存在するコンピュータプログラムの製品を提供する。
【発明を実施するための最良の形態】
【００３９】
以下、非限定的な例を用いて、添付の図面を参照しながら、本発明を具体化するシステムについて説明する。
【００４０】
図１−図５に関連するスレンダボディフレキシブルシステムの記載は、具体的には、制限された（limited）自由度の問題に対処するものではなく、従来技術インクリメンタルトータルラグランジュ法（incremental total Lagrangian formulation）における共回転アップデーティング（co-rotational updating）に対処するものであるが、これは、ソルバの一般的な操作を理解するのに役立つ。明確にするために、ＮＡＦＥＭＳ  Ｗｏｒｌｄ  Ｃｏｎｇｒｅｓｓの論文（以後、ＮＡＦＥＭＳの論文とする）に開示されたプロセスについて、ここで再び説明する。本記載の第１の部分として、両端におけるフレキシブルシステムの予め定められた基準変位が、完全に決定される場合について検討する。
【００４１】
このプロセスの目的は、予め定められた変位境界条件を受けたときの、フレキシブルシステムの挙動を判定することである。図１における参照番号２は、フレキシブルシステムの開始コンフィギュレーションを示している。システムの両端には、リファレンスフレーム４および６が示されている。これらのリファレンスフレームは、ローカルリファレンスフレームであり、両端におけるフレキシブルシステムの向きを表している。したがって、単語「変位」は、フレキシブルシステムの端の位置だけでなく、フレキシブルシステムの向きも含むことを理解されたい。したがって、点の変位には、６つの自由度が含まれる。
【００４２】
図１における参照番号１２は、フレキシブルシステムの最終コンフィギュレーションを示しているのに対し、参照番号１４および１６は、フレキシブルシステムの両端におけるリファレンスフレームを示している。第１の端（リファレンスフレーム６および１６）では、変位は、位置の変位を含み、リファレンスフレーム６および１６の軸は、平行である。第２の端（リファレンスフレーム４および１４）では、変位は、位置および向きの変位を含み、リファレンスフレーム１４の軸は、リファレンスフレーム４の軸と比較すると、向きが変わっている。
【００４３】
図１は、フレキシブルシステムの中間コンフィギュレーション２０、２２および２４をさらに示している。図１における中間コンフィギュレーションの数は、説明の目的上選択されたに過ぎない。一コンフィギュレーションから次のコンフィギュレーションに変化させる（これが１回のインクリメンタルステップを表す。）ために、予め定められた変位の一部がフレキシブルシステムに与えられる。３つの中間コンフィギュレーションを伴う、したがって、４回のインクリメンタルステップを伴う図１の例では、フレキシブルシステムに予め定められた変位の４分の１の変位を与えることによって、一コンフィギュレーションから次のコンフィギュレーションに変化する。例えば、第１インクリメンタルステップにおいて、フレキシブルシステムの端に、予め定められた変位の４分の１の並進変位および回転変位を与えることによって、開始コンフィギュレーション２から第１中間コンフィギュレーション２０に変化する。
【００４４】
フレキシブルシステムの各中間コンフィギュレーションおよび最終コンフィギュレーションは、幾何学的非線形有限要素解析を用いて計算される。具体的には、フレキシブルシステムは、一連の有限要素として表現される。説明の目的上、第１インクリメンタルステップについて検討する。この第１インクリメンタルステップにおいて、開始コンフィギュレーション２から第１中間コンフィギュレーション２０に変化する。図２は、リファレンスフレーム４と、いくつかの要素とを伴う開始コンフィギュレーション２の部分図を示している。これは、複数の要素３０、３５、４０、４５を示している。各要素は、節点３１、３６、４１、４６および５１の２つの節点の間を伸びる。節点３１は、実際には、フレキシブルシステムの端である。各要素３０、３５、４０、４５は、図２において３２、３７、４２および４７の参照番号が付されたエレメンタリリファレンスフレームと関連付けられる。図２は、第１中間コンフィギュレーション２０および各節点の部分図をさらに示している。要素およびリファレンスフレームには、３０だけ増加した参照番号が付されている。
【００４５】
各インクリメンタルステップにおいて、予め定められた変位の一部を受けたときのフレキシブルシステムのトータルポテンシャルエネルギを最小化する。フレキシブルシステムのトータルポテンシャルエネルギは、フレキシブルシステムを形成する各有限要素のトータルポテンシャルエネルギの合計として計算される。これは、ＮＡＦＥＭＳの論文に開示された反復プロセスによって実行される。この反復プロセスはまた、以下で説明するように実行することができる。
【００４６】
各インクリメンタルステップにおいて、共回転アップデートを伴うトータルラグランジュフォーマリズム（total Lagrangian formalism）を用いる。ＮＡＦＥＭＳの論文のセクション４．３．２．１にて論じられているように、所与のインクリメンタルステップ（または、ＮＡＦＥＭＳの論文にて用いられている表現「palier」）中の変位は、前のインクリメンタルステップにおいて得られた解に対応するコンフィギュレーションに関連して、測定される。すなわち、一インクリメンタルステップにおいて、フレキシブルシステムのトータルポテンシャルエネルギの最小化は、インクリメンタルステップ開始時のフレキシブルシステムのコンフィギュレーションにおける要素に関連付けられたエレメンタリリファレンスフレームを用いて、各要素ごとに、実行される。図２の例において、インクリメンタルステップ開始時のコンフィギュレーションは、開始コンフィギュレーション２である。したがって、第１中間コンフィギュレーションを得るために、開始コンフィギュレーションに対応するリファレンスフレーム、すなわち図２のエレメンタリリファレンスフレーム３２、３７、４２および４７を用いて、反復的にフレキシブルシステムのトータルポテンシャルエネルギを最小化する。
【００４７】
しかしながら、次のインクリメンタルステップでは、第１中間コンフィギュレーションに対応する最新のリファレンスフレームを考慮する。すなわち、次のインクリメンタルステップは、第１中間コンフィギュレーション２０から第２中間コンフィギュレーション２２に変化させる。第２インクリメンタルステップの計算は、第１中間コンフィギュレーションに関連付けられたエレメンタリリファレンスフレーム６２、６７、７２および７７を用いて実行される。これにより、次のインクリメンタルステップの前に、各インクリメンタルステップの最後に新たなリファレンスフレームを計算することが必要となる。ＮＡＦＥＭＳの論文には、各要素ごとの最新のリファレンスフレームの計算についての詳細が提供されていない。すなわち、共回転アップデートが提案されてはいるが、この共回転アップデートの実行に関しては、説明がなされていない。
【００４８】
本発明は、所与の一有限要素に関して、有限要素の節点における計算された変位ベクトルのコンポーネントに基づいて、変位ベクトルの曲げ回転コンポーネントを考慮することなく、共回転アップデートを実行することを提案する。
【００４９】
本発明の解決策は、最新のエレメンタリリファレンスフレームにおいて、軸の１つが有限要素の節点を結ぶ直線と平行であることを保証するという利点を有する。この特性は、各インクリメンタルステップにおいて維持され、この特性により、すべての計算が類似したものとなる。したがって、それらの計算は、より単純、かつ、より容易になる。
【００５０】
図３は、図２の有限要素４０および７０の拡大図である。これは、開始コンフィギュレーションにおける節点４１および４６を示しており、有限要素４０は、節点４１と節点４６との間を伸びる。有限要素４０のエレメンタリリファレンスフレーム４２は、軸ｘ、ｙおよびｚを有する。ｘ軸は、節点４１と節点４６とを結ぶ直線に平行である。
【００５１】
フレキシブルシステムのトータルポテンシャルエネルギを最小化することによって、フレキシブルシステムの各節点ごとの変位ベクトルが提供される。変位ベクトルは、並進コンポーネント（ｕｘ，ｕｙ，ｕｚ）および回転コンポーネント（θｘ，θｙ，θｚ）から構成される。並進コンポーネント（ｕｘ，ｕｙ，ｕｚ）は、各節点の位置の変位を表している。並進コンポーネントは、容易に理解されるものであり、図３の例では、節点４１から節点７１への位置の変位、または、節点４６から節点７６への位置の変位に相当する。
【００５２】
回転コンポーネント（θｘ，θｙ，θｚ）は、各節点の向きの変位を表している。この回転コンポーネントは、システムがフレキシブルシステムであるため、所与の節点におけるシステムの向きが一インクリメンタルステップごとに変化することを表している。
【００５３】
回転コンポーネントは、曲げ回転コンポーネントを含む。曲げ回転コンポーネントは、インクリメンタルステップ終了時に得られたコンフィギュレーション内の所与の節点におけるフレキシブルシステムに与えた曲げ（bending）を表している。
【００５４】
回転コンポーネントは、ねじり回転コンポーネントをさらに含む。ねじり回転コンポーネントは、所与の節点におけるフレキシブルシステムに与えたねじり（torsion）を表している。
【００５５】
節点４１における回転コンポーネント（θｘ，θｙ，θｚ）が、要素４０のリファレンスフレーム４２を用いて表現される場合（ただし、ｘ軸は、節点４１から節点４６へ伸びる直線に平行である）、ねじりコンポーネントは、ｘ軸に沿って測定されたコンポーネントθｘである。曲げ回転コンポーネントは、コンポーネントθｙおよびθｚである。すなわち、曲げ回転コンポーネントは、エレメンタリリファレンスフレームのｘ軸に直交する軸に沿った節点の回転を表していて、エレメンタリリファレンスフレームのｘ軸は、有限要素の節点を結ぶ直線に平行である。
【００５６】
ねじり回転コンポーネントはθｘであり、これは、エレメンタリリファレンスフレームのｘ軸に沿って節点において与えたねじりを表している。
【００５７】
第１の実施例において、最新のエレメンタリリファレンスフレームは、有限要素の節点における変位ベクトルの並進コンポーネントのみに基づいて計算される。この実施例は、２次元でのみ変形するフレキシブルシステムに適する。このような場合、フレキシブルシステムは、ねじり回転制約を受けず、２Ｄ平面内で曲がるだけである。このような平面変形（plane deformation）の例は、両端が固定され、重力の作用を受け、かつ、一端において強制的な水平の圧縮変位を受ける水平梁によって提供される。
【００５８】
この場合、次のように、最新のエレメンタリリファレンスフレームを計算することができる。インクリメンタルステップ開始時のコンフィギュレーションにおいて、２つの節点Ｐ１とＰ２とを結ぶ有限要素Ｐ１Ｐ２について検討する。この有限要素用のエレメンタリフレームは、３つの軸ｘ、ｙおよびｚを有する。ｘ軸は、２つの節点を結ぶ直線（Ｐ１Ｐ２）に平行である。ｙ軸は、２Ｄ平面内でｘ軸に直交する。ｚ軸は、２Ｄ平面に直交する。
【００５９】
（ｕｘ１，ｕｙ１，ｕｚ１）を、有限要素のエレメンタリリファレンスフレームを用いて測定された、有限要素の第１の節点Ｐ１における計算された変位ベクトルの並進コンポーネントとする。（θｘ１，θｙ１，θｚ１）を、有限要素の第１の節点Ｐ１における計算された変位ベクトルの回転コンポーネントとする。同様に、（ｕｘ２，ｕｙ２，ｕｚ２）および（θｘ２，θｙ２，θｚ２）を、有限要素のエレメンタリリファレンスフレームを用いて測定された、第２の節点Ｐ２に関する変位ベクトルの並進コンポーネントおよび回転コンポーネントとする。Ｐ１’およびＰ２’は、インクリメンタルステップにおいて計算された、変位した節点である。
【００６０】
２次元の場合のフレキシブルシステムは、ねじり制約を有さないので、θｘ１＝θｘ２＝０である。第１の節点および第２の節点におけるねじり回転制約はゼロである。フレキシブルシステムは２Ｄ平面内を動くので、ｕｚ１＝ｕｚ２＝０である。
【００６１】
したがって、第１の節点および第２の節点における曲げ回転コンポーネントは、θｚ１およびθｚ２によって表現される。
【００６２】
この実施例において、最新のエレメンタリリファレンスフレームは、次のように計算される。エレメンタリリファレンスフレームのｘ軸は、変位した節点Ｐ１’とＰ２’とを結ぶ直線（Ｐ１’Ｐ２’）に平行であるように設定される。ｘ軸は、Ｐ１’からＰ２’へ伸びる方向に向けられる。ｙ軸は、フレキシブルシステムが動いている２Ｄ平面内でｘ軸に直交する。ｚ軸は、この場合もまた、２Ｄ平面に直交する。最新のエレメンタリリファレンスフレームの計算に関して、変位ベクトルの回転コンポーネントは考慮されない。
【００６３】
第２の実施例において、最新のエレメンタリリファレンスフレームは、有限要素の節点における変位ベクトルの並進コンポーネント、および、ねじり回転コンポーネントに基づいて計算される。この実施例は、３次元で変形するフレキシブルシステムに適する。これは、例えば、ケーブルに与えたねじり制約がケーブルの３Ｄ変形をもたらすケーブルに適用され、より一般的には、直線状のニュートラルラインが非平面曲線（non-planar curve）に変形する任意の状況に適用される。
【００６４】
節点Ｐ１およびＰ２に関する変位ベクトルは、第１の実施例と同一の記号が用いられる。この実施例では、θｘ１＝θｘ２＝０、または、ｕｚ１＝ｕｚ２＝０と仮定する理由はない。最新のエレメンタリリファレンスフレームは、次のように計算される。まず、節点Ｐ１およびＰ２のねじりコンポーネントθｘ１およびθｘ２について検討する。最新のエレメンタリリファレンスフレームは、以下によって計算される。
−ｘ軸に沿って（すなわち、節点Ｐ１とＰ２とを結ぶ直線（Ｐ１Ｐ２）に沿って）、エレメンタリリファレンスフレームに、節点Ｐ１およびＰ２のねじり回転コンポーネントの平均（θｘ１＋θｘ２）／２に等しい角度分の回転を与えること、
−回転したフレームに、節点Ｐ１から節点Ｐ２に向かうユニタリベクトル（unitary vector）を、節点Ｐ１’から節点Ｐ２’に向かうユニタリベクトルに変える最小回転に等しい回転を与えること。この回転は、変位ベクトルの並進コンポーネント（ｕｘ１，ｕｙ１，ｕｚ１）および（ｕｘ２，ｕｙ２，ｕｚ２）に基づいて計算され、変位ベクトルの回転コンポーネント（θｘ１，θｙ１，θｚ１）および（θｘ２，θｙ２，θｚ２）を考慮に入れない。
【００６５】
すなわち、まず、インクリメンタルステップ中に使用されるエレメンタリリファレンスフレームが、節点に与えたねじり回転コンポーネントの平均に一致する量だけ回転される。次いで、最新のエレメンタリリファレンスフレームのｘ軸が、変位した節点Ｐ１’とＰ２’とを結ぶ直線（Ｐ１’Ｐ２’）に平行となるように、回転したフレームがさらに回転される。
【００６６】
第１の実施例のように、最新のリファレンスフレームのｘ軸は、変位した節点Ｐ１’とＰ２’とを結ぶ直線に平行である。しかしながら、第１の実施例とは異なり、節点に与えたねじり変位を考慮に入れるため、ｙ軸およびｚ軸は回転されるので、ｘｙ平面が一定に保たれるという理由はない。
【００６７】
双方の実施例において、最新のエレメンタリリファレンスフレームは、有限要素の節点における変位ベクトルの曲げ回転コンポーネントを考慮することなく、計算される。また、双方の実施例において、エレメンタリリファレンスフレームのｘ軸は、有限要素の節点を結ぶ直線に平行なままである。
【００６８】
上記にて開示したプロセスは、各有限要素ごとに、最新のリファレンスフレームを提供する。この最新のリファレンスフレームは、次のインクリメンタルステップにおける反復計算のために使用することができる。エレメンタリリファレンスフレームのアップデートは、有限要素解析において用いられる近似が終始妥当であることを保証する。具体的には、
【００６９】
【数４】
2006294038.tif 000006
【００７０】
となるように、回転θｙおよびθｚが十分に小さい間は、ＮＡＦＥＭＳの論文のセクション４．１．１．１に開示された、ケーブルキネマティックス（cable kinematics）  
ｕ（ｘ，ｙ，ｚ）＝ｕ０（ｘ）−ｙ．θｚ（ｘ）＋ｚ．θｙ（ｘ）
ｖ（ｘ，ｙ，ｚ）＝ｖ０（ｘ）−ｚ．θｙ（ｘ）
ｗ（ｘ，ｙ，ｚ）＝ｗ０（ｘ）＋ｙ．θｚ（ｘ）
は、終始妥当である。エレメンタリリファレンスフレームが上述のようにアップデートされるとき、この仮定は妥当である。
【００７１】
エレメンタリリファレンスフレームの共回転アップデートは、各節点における変位ベクトルの曲げコンポーネントを考慮に入れない。変位ベクトルの曲げコンポーネントは、節点リファレンスフレームをアップデートするプロセスによって、考慮に入れられる。
【００７２】
加えて、インクリメンタルステップは、フレキシブルシステムの有限要素表現を提供する。これは、システムの傾斜不連続な（slope-discontinuous）表現である。フレキシブルシステムのスムーズな（、かつ傾斜不連続な）表現を得るために、本発明は、節点リファレンスフレームを使用することを提案する。節点リファレンスフレームは、フレキシブルシステムの各節点に関連付けられる。節点フレームは、変位ベクトルの曲げ回転コンポーネントを考慮に入れることによって、エレメンタリリファレンスフレームとは別にアップデートされる。
【００７３】
図３は、節点リファレンスフレームをさらに示している。開始コンフィギュレーション２において、節点リファレンスフレームは、以下のように定義される。
−各節点リファレンスフレームのｘ軸は、フレキシブルシステムに接する。
−各節点リファレンスフレームのｙ軸は、ニュートラルラインに沿ったケーブルシステムのねじりを表すものである。すなわち、フレキシブルシステムの一端における、例えば、第１の節点３１における、任意の所与の方向のｙ軸が選択される。次の節点に関して、ｙ軸は、ｙ軸に直交するものとして定められ、フレキシブルシステムに与えたねじり制約を考慮に入れる。したがって、ある平面において、フレキシブルシステムがねじりを受けない場合、様々な節点リファレンスフレーム用のｙ軸は、すべてその平面内に位置するか、または、すべてその平面に直交するかのいずれかである。
−各節点リファレンスフレームのｚ軸は、フレームのすべての軸が直交するような軸とする。
【００７４】
各インクリメンタルステップにおいて、節点リファレンスフレームは、次のようにアップデートされる。節点Ｐ１における節点リファレンスフレームは、節点における変位ベクトルの回転コンポーネント（θｘ１，θｙ１，θｚ１）分だけ回転されることによって、アップデートされる。
【００７５】
最新の節点リファレンスフレームの基点（origin）は、変位した節点Ｐ１’である。これは図３に示されており、節点リファレンスフレーム４３および４８が、それぞれ、節点４１および４６に対して用いられ、さらに、節点リファレンスフレーム７３および７８が、節点７１および７６に対して用いられる。図３はさらに、２点破線によって、フレキシブルシステムのスムーズコンフィギュレーション（smoothed configuration）を示している。節点リファレンスフレームのアップデートは、変位ベクトルの曲げ回転コンポーネント、および、ねじり回転コンポーネントを考慮に入れる。
【００７６】
節点リファレンスフレームが提供されることによって、各インクリメンタルステップにおいて、フレキシブルシステムの形状を計算することが可能となる。フレキシブルシステムの形状は、節点の位置と、節点リファレンスフレームとに基づいて計算される。具体的には、フレキシブルシステムの形状は、節点の位置と、様々な節点における節点リファレンスフレームとを用いて、補間される。以下の補間のルールを用いることができる。
−フレキシブルシステムは、各節点を通らなければならない。
−各節点において、フレキシブルシステムは、節点リファレンスフレームのｘ軸に接しなければならない。
−各節点において、ｙ軸およびｚ軸は、ケーブルに与えたねじり制約を表すものである。
【００７７】
フレキシブルシステムの形状を補間するために、３次元スプライン補間などの、最先端の補間ツールを例えば使用することができる。
【００７８】
補間するステップは、好ましくは、フレキシブルシステムの変形した形状をユーザに表示する前に実行される。したがって、ユーザには、フレキシブルシステムの有限要素ビュー（finite element view）ではなく、変形したフレキシブルシステムのスムーズビュー（smoothed view）を提供することができる。
【００７９】
各インクリメンタルステップにおいて、節点リファレンスフレームによって、ステップ開始時にフレキシブルシステムに与えた初期応力が考慮され得る。上記にて説明したように、エレメンタリリファレンスフレームは、有限要素の節点に与えた変位ベクトルの曲げ回転コンポーネントを考慮することなく、アップデートされる。それに対し、節点リファレンスフレームは、節点に与えた変位ベクトルのすべての回転コンポーネントに基づいて、アップデートされる。
【００８０】
エレメンタリリファレンスフレームおよび節点リファレンスフレームをアップデートした後、一節点における、（その節点の他の一辺上の）エレメンタリリファレンスフレームと、節点リファレンスフレームとの間の角度差は、曲げ回転コンポーネント量を表している。この差はまた、フレキシブルシステムに与えた初期応力量を表している。
【００８１】
本発明はまた、インクリメンタルステップ開始時にフレキシブルシステムに与えた初期応力を測定する手段として、節点リファレンスフレームを使用することを提案する。これは、フレキシブルシステムの各有限要素の各節点において、節点リファレンスフレームと、有限要素のエレメンタリリファレンスフレームとを比較することによって、実行される。この比較（すなわち、節点リファレンスフレームと、エレメンタリリファレンスフレームとの差）は、フレキシブルシステムに与えた初期応力を表している。１つの可能な比較は、エレメンタリリファレンスフレームを節点リファレンスフレームに変換する回転を計算することにある。この回転の角度は、節点の初期応力を表していて、インクリメンタルステップ中に、最小化されるエネルギに対する、対応するコントリビューション（contribution）を計算するために、使用することができる。所与の一節点における、その節点を含む２つの有限要素用のエレメンタリリファレンスフレームと、節点リファレンスフレームとの差を考慮してもよいし、あるいは、所与の要素用のエレメンタリリファレンスフレームと、その要素を形成する節点における節点リファレンスフレームとの差を考慮してもよい。差がフレキシブルシステム全体にわたって足し合わされたとき、双方の例は、同じ結果になるか、または、比例する結果となる。
【００８２】
例えば、ｘｙ平面内の平面運動については、エレメンタリひずみエネルギ（elementary strain energy）に関する以下の式を用いることができる。
【００８３】
【数５】
2006294038.tif 000007
【００８４】
この式において、上付き文字が付された項（superscripted term）は、節点リファレンスフレームと、インクリメンタルステップ開始時の対応するエレメンタリリファレンスフレームとの間の、要素節点における角度差を表している。上付き文字が付されていない項（non-superscripted term）は、現行の反復での節点変位を表している。
【００８５】
最小化されるエネルギについて、節点リファレンスフレームと、エレメンタリリファレンスフレームとの差を表す項を使用することにより、フレキシブルシステムの挙動を計算する際に、不安定性が生じる場合がある。このような不安定性は、一インクリメンタルステップから次のインクリメンタルステップに進んだ際の、フレキシブルシステムの形状の大幅な変化として生じる。さらに、このような項を用いることによって、インクリメンタルステップ中に収束解を見つけることが、より難しくなる場合がある。
【００８６】
こうした問題を回避または制限するため、対応する１つまたは複数の項に、スケーリング係数μＩＳを乗じることができる。スケーリング係数は、０以上１未満である。エレメンタリひずみエネルギに関する例示的な式は以下のようになる。
【００８７】
【数６】
2006294038.tif 000008
【００８８】
スケーリング係数μＩＳが１であるとき、初期応力は考慮される。スケーリング係数μＩＳが０であるとき、初期応力は考慮されない。スケーリング係数μＩＳは、好ましくは、結果に応じて、インクリメンタルステップ全体を通して、変更される。インクリメンタルステップ中に収束解を見つけることができないとき、または、フレキシブルシステムが不安定であると分かったとき、通常スケーリング係数を小さくするであろう。ユーザがこれを実行することができる。代替として、例えば、インクリメンタルステップにおける現行の反復回数に応じて、スケーリング係数を自動的に適合させることができる。別の実施例においては、スケーリング係数は、要素の平均節点リファレンスフレーム（mean nodal frame of reference）（要素の端節点における節点リファレンスフレームの平均）と、リファレンスエレメンタリフレームとの差の値によって決まる。この実施例では、微小な差がインクリメンタルステップにおける反復プロセスの結果の収束も示すという事実が用いられる。すなわち、差が小さいほど、有限要素はフレキシブルシステムの挙動をより良くシミュレートする。
【００８９】
図４および図５は、インクリメンタルステップ終了時の有限要素の例を示している。説明の目的上、いずれの図も２次元であり、これらは、有限要素７０と、節点フレーム７３および７８と、エレメンタリフレーム７２とを示し、さらに、平均節点フレーム８０も示している。上記にて説明したように、平均節点フレームは、節点フレームの平均として計算される。その基点は、節点間の中央であり、したがって、エレメンタリフレームと同一の基点である。平均節点フレームの各ベクトルは、節点フレームの対応するベクトルの正規化した和である。図４において、エレメンタリフレーム７２と、平均節点フレーム８０との間の角度差は小さい（この例では、数度である）。これは、良好な収束を表し、スケーリング係数μＩＳの適切な選択を表している。図５において、エレメンタリフレーム７２と、平均節点フレーム８０との間の角度差は重大であり（約４５°）、これは、スケーリング係数μＩＳの値が大きすぎることを表している。３次元の例では、平均節点フレームと、エレメンタリリファレンスフレームとの間の差は、エレメンタリリファレンスフレームを平均節点フレームに変換する回転の角度として計算することができる。この角度の値が１０°未満であれば、スケーリング係数が適切な値であることを表し、適切な収束であることを表す。
【００９０】
上述したように、各インクリメンタルステップにおいて、平衡解を求めるために反復サーチが実行される。各反復は、２段階のプロセスからなる。このプロセスの第１段階は、エネルギ表面に対する降下方向を見つけることにある。このプロセスの第２段階は、最適ステップを見つけるために、降下方向をラインサーチすることにある。最適ステップを見つけるには、収束判定条件（stopping criterion）が必要である。ＮＡＦＥＭＳの論文のセクション４．３．２．４では、各反復において、余剰ベクトル（residual vector）およびトータルエネルギを求めることが提案されている。次いで、個々の収束判定が、余剰ベクトルのノルム、および、トータルエネルギ変動のモノトニ（monotony）に対して行われる。余剰ベクトルのノルム、および、エネルギ変動のモノトニの両方に基づく、この収束判定条件は、有効である。しかしながら、それほど厳しくはない収束判定条件を与えることが有用な場合もある。これは、平衡を見つける前に、フレキシブルシステムが不安定状態を抜け出す場合、特に有用である。
【００９１】
本発明は、余剰変動と変位変動とのスカラ積として定義される余剰エネルギノルムを、収束判定条件として使用することを提案する。このような収束判定条件を使用することによって、システムは、より容易に不安定状態を抜け出すことができる。システムが分岐を経て、所与の一コンフィギュレーションから非常に異なるコンフィギュレーションに否応なく変化するときに、これは特に明らかである。提案した収束判定条件を使用することによって、各インクリメンタルステップにおいて、エネルギの最小化を反復して実行する間に、解を見つけることがより容易になる。したがって、計算は、より短くなる。この収束判定条件の使用は、上記にて開示した共回転アップデートには依拠しない。
【００９２】
したがって、これまでのプロセスの説明において、両端におけるフレキシブルシステムの予め定められた基準変位は、完全に決定される。すなわち、フレキシブルシステムの位置および向きは、完全に制約される。本発明では、フレキシブルシステムが、その端においても、いくつかの自由度を有する場合について検討する。フレキシブルシステムに関する一例は、フレキシブルシステムの一点（しばしば、その一端）が、位置については固定されているが、回転については自由である場合である。これは、その点が、並進については固定されているが、回転については３自由度を有することを意味する。別の例は、フレキシブルシステムに関するすべり点（sliding point）である。システムは、一方向に沿って自由にすべることができ、おそらくはすべり方向のまわりを自由に回転することができる。この例では、フレキシブルシステムは、並進について１自由度を有し、回転について１自由度を有する。
【００９３】
この問題は、非線形の構造システムであるスレンダボディフレキシブルシステムに関して存在する。この問題はまた、比較的硬い、もしくは、非フレキシブルな梁またはシェルなどの線形の構造システムに関しても存在する。
【００９４】
この問題はまた、線形の連続システムまたは非線形の連続システムに関しても存在する。線形の連続システムには、例えば、巨大な飛行機の固定具、エンジンキャップ、シリンダブロック、クランクシャフト、または、複雑な形状の巨大な機械部品などを通常のＣＡＤを用いてモデル化したものなどの複雑なソリッド形状が含まれる。非線形の連続システムには、自動車産業において使用されるゴム製のサイレントブロック、または、ウィンドサーフィンのマスト用に使用されるジョイント材などの巨大な変形しやすいオブジェクトが含まれる。各有限要素に関する変位ベクトルは、並進コンポーネントしか有さないが、巨視的オブジェクトのモデル化された変形した形状にはおそらく、概して、回転が含まれるであろう。
【００９５】
線形システムに関しては、従来技術において、解法が存在する。このような解法は、「Analyse des structures par elements finis, J.F. Imbert, CEPADUES Editions, 1979-1984, page 276」、または、「MSC NASTRAN Version 64, Handbook for linear Analysis, 1985, The MacNeal-Schwendel Corp., page 4.4.1, Chapter 4.4, Multipoint Constraint Operations」に開示されている。これらの文献は、制約をファクトライズすること、すなわち、制約または自由度の一部を他の自由度の関数として表現することを提案している。具体的には、Ｘを、所与の線形システムにおける、自由度の列ベクトルとする。解くべき一般的問題は、以下のように記述することができる。    
Ｋ．Ｘ＋ＡＴ．λ＝Ｆａｐｐ  かつ  Ａ．Ｘ＝ｂ          （１）
ただし、Ｋは剛性行列であり、Ａは変位制約のヤコビ行列であり、λはラグランジュ乗数の列ベクトルであり、ｂは右側制約（constraint right-hand sides）の列ベクトルであり、Ｆａｐｐはシステムに与えた力の列ベクトルである（ＡＴは行列Ａの転置を示す）。従来技術による解法は、制約をファクトライズすること、すなわち、
Ｘ＝Ｎ．Ｘａ．＋Ｍ          （２）
となるような行列ＮおよびＭとともに、自由度の独立サブセット（independent subset）Ｘａを見つけることにある。
【００９６】
これにより、式（１）におけるＸの値を置き換えることが可能になる。Ａ．Ｎ＝０と仮定すると、ＮＴを左から乗じたときの式（１）により、
ＮＴ．Ｋ．Ｎ．Ｘａ＝ＮＴ．（Ｆａｐｐ−Ｋ．Ｍ）
が与えられる。従来技術の解法を用いることによって、この式を独立自由度Ｘａについて解くことができる。Ｘａ（これは、緩和された自由度（relaxed degree of freedom）のセットを含む）が見つかると、式（２）によって、完全な変位を復元することができる。
【００９７】
自由度のセットを独立サブセットおよび依存サブセット（dependent subset）に自動的に分割することは、線形問題の場合でさえも時間を要するプロセスであるということに加えて、非線形問題の場合には、この解法によって、いくつかのさらなる問題が生じる。第１に、非線形システムでは、上記にて説明したように、リファレンスフレームはアップデートされる必要がある。これは、リファレンスフレームがアップデートされるたびに、制約が繰り返しファクトライズされなければならないことを意味する。すなわち、各インクリメンタルステップ後に、式（２）における行列Ｎの値は、再度計算される必要がある。第２に、数字が限られた桁数で記憶されるコンピュータシステムでは、制約をファクトライズするのに必要とされる逆行列化によって、計算のロバスト性（robustness）が低下する場合がある。最後に、剛性行列Ｋは、各インクリメンタルステップごとに変わる可能性があり、これによって、ＮＴ．Ｋ．Ｎを繰り返し計算することが必要とされる。これは、時間を要するものであることが分かる。
【００９８】
線形システムでは、共回転アップデートは必要とされない。それでもやはり、本発明は、制約をファクトライズする必要性を回避し、単純なリファレンスフレームの変換のみを必要とする。
【００９９】
したがって、リファレンスフレームのアップデートを伴う非線形システムにおいて、制約が緩和されることを可能にするための要求が存在する。本発明は、部分的に緩和された自由度と、部分的にインポーズされた（imposed）自由度とからなる制約を伴うシステムの節点用のローカルリファレンスフレームを用いてエネルギを最小化するために、ハイブリッド（hybrid）リファレンスフレームを使用することを提案する。
【０１００】
解放された自由度がローカルリファレンスフレームの軸に沿って表現されるように、ローカルリファレンスフレームは選択される。次いで、式（１）において、緩和された（インポーズされていない（non-imposed））自由度に対応する行（line）を単に選択することによって（換言すれば、インポーズされた自由度、すなわち、緩和されていない自由度に対応する行および列を剛性行列から取り除くことによって）、制約されていない自由度のサブセットＸａを見つけることができる。
【０１０１】
これは、上述した問題に対する解決策を提供する。本発明の解決策は、各インクリメンタルステップにおいて、特に選択されたリファレンスフレームの変換を行うことのみを必要とする。しかしながら、いずれにせよ、リファレンスフレームの変換は、エネルギを最小化するために、各ステップにおいて必要である。以下で説明するように、ハイブリッドリファレンスフレームのセットへの変換は、計算の観点からみると、それほど複雑ではない。
【０１０２】
図６を参照して、この解決策を例示する。図６は、２つの有限要素３５および４０を示しており、これらは、それぞれ、節点３６、４１、および、４１、４６を伴う。図６は、第１の有限要素３５および第２の有限要素４０用のエレメンタリリファレンスフレーム３７および４２をさらに示している。単純化するために、フレーム３７および４２は、それぞれ、Ｒ１およびＲ２と表記されている。Ｘ１を点３６における自由度の列ベクトルとするのに対し、Ｘ２およびＸ３は、それぞれ、点４１および４６における自由度の列ベクトルである。
【０１０３】
緩和された自由度が存在しない場合、まず、剛性行列が、各有限要素３５および４０用に、エレメンタリリファレンスフレーム、または、有限要素に関連する任意のフレームを用いて計算される。剛性行列が、各要素用に、対応するエレメンタリフレームを用いて計算されると、様々な行列が、グローバルリファレンスフレームを用いて表現され、次いで、以下の式（３）ないし（６）において説明するように、これらを足し合わせることができる。
【０１０４】
具体的には、第１の要素３５に関して、ローカルフレームＲ１を用いて、平衡方程式Ｋ．Ｘ＝Ｆが記述される。
【０１０５】
【数７】
2006294038.tif 000009
【０１０６】
同様に、第２の有限要素に関して、第２の要素用の対応するローカルフレームＲ２を用いて、
【０１０７】
【数８】
2006294038.tif 000010
【０１０８】
が得られる。ただし、各Ｘｋが６自由度を有する場合、各Ｋｉｊは、実際には６×６の行列である。式（３）および（４）において、行列に対する添え字（subscript）１および２は、行列がローカルリファレンスフレームＲ１およびＲ２を用いて表現されているということを表している。
【０１０９】
制約がまったく緩和されていない場合、図６において８２の参照番号が付されたグローバルリファレンスフレームＲＧを用いてこれらがすべて表されるように、リファレンスフレームの変換が式（３）および（４）に単に適用される。ローカルリファレンスフレームＲ１をグローバルリファレンスフレームＲＧに変換するための行列が、Ｐ１と表記されるとすると、
【０１１０】
【数９】
2006294038.tif 000011
【０１１１】
が得られる。ただし、添え字Ｇは、自由度がグローバルリファレンスフレームＲＧを用いて表現されていることを示す。したがって、式（３）は、
【０１１２】
【数１０】
2006294038.tif 000012
【０１１３】
に変形される。ただし、Ｐ１は、ローカルリファレンスフレームＲ１をグローバルリファレンスフレームＲＧに変換するための、リファレンスフレーム変換用の行列である。したがって、式（３）は、フレームＲＧを用いて、
【０１１４】
【数１１】
2006294038.tif 000013
【０１１５】
と表されるのに対し、式（４）は、フレームＲＧを用いて、
【０１１６】
【数１２】
2006294038.tif 000014
【０１１７】
と表される。Ｐ２は、ローカルリファレンスフレームＲ２をグローバルリファレンスフレームＲＧに変換するための、リファレンスフレーム変換用の行列である。すべての有限要素について式（５）および（６）を足し合わせることによって、グローバルな平衡方程式が得られる。２つの有限要素のみを伴う図６の例では、これにより、以下の式が与えられる。
【０１１８】
【数１３】
2006294038.tif 000015
【０１１９】
これらの式には、いくらかのあいまいさが存在する。例えば、Ｆ２は、Ｒ１を用いて表現された、第２の節点４１に与えた力を示すと同時に、Ｒ２を用いて表現された、第２の節点４１に与えた力を示す。このあいまいさは、行列Ｐ１およびＰ２の使用の結果、生じる。したがって、単純化の目的上、あいまいさを回避するためのさらなる添え字は使用しない。式（６）は、グローバルリファレンスフレームを用いて計算された、剛性行列の形式を示している。図６の例は、各ローカルリファレンスフレームを用いて計算された剛性行列を足し合わせることによって、３以上の有限要素に対して一般化される。必要とされるリファレンスフレームの共回転アップデートにより、フレキシブルシステムの各要素ごとのフレーム変換Ｐｉの計算が、各インクリメンタルステップにおいて実行される。
【０１２０】
本発明は、エネルギを最小化するために、ハイブリッドリファレンスフレームを使用することを提案する。いくつかのインポーズされた自由度と、いくつかの解放された自由度とを有する節点に対し、グローバルフレームを使用する代わりに、ローカルフレームを使用する。自由度がローカルフレームの軸に沿って表現されるように、ローカルフレームは選択される。
【０１２１】
例えば、節点４１において、少なくとも１つの自由度が解放されていると仮定されたい。そうでなければ節点４１は固定される。ＲＬを、その解放された自由度が１つの軸に沿って表現されるリファレンスフレームとする。このフレームは、図６では、８３の参照番号が付されている。例えば、緩和された制約または自由度が、タンジェンシャルな（tangential）ものであり、それによってシステムが一方向に自由にすべることができる場合、タンジェンシャルな節点フレームを、ローカルリファレンスフレームＲＬとして使用することができる。システムが一方向に沿って自由にすべることができ、かつ、同一方向のまわりを自由に回転することができるならば、同一の節点フレームを使用することができる。このような場合、行列Ｘ１に関しては、同一のグローバルリファレンスフレームＲＧに変換する、フレーム変換Ｐ１を使用する。同様に、行列Ｘ３に関しては、同一のグローバルリファレンスフレームＲＧに変換する、フレーム変換Ｐ２を使用する。しかしながら、行列Ｘ２に関しては、フレームはローカルフレームＲＬに変換される。上記の式（３）は、
【０１２２】
【数１４】
2006294038.tif 000016
【０１２３】
に変形される。ただし、最後の添え字は、該当する場合、Ｘ１およびＸ２を表現するために使用されるリファレンスフレームを示し、
−Ｐ１Ｌは、リファレンスフレームを、ローカルリファレンスフレームＲ１からローカルリファレンスフレームＲＬに変換するための行列を示し、
−Ｐ２Ｌは、リファレンスフレームを、ローカルリファレンスフレームＲ２からローカルリファレンスフレームＲＬに変換するための行列を示す。
ローカルフレームＲＬがローカルフレームＲ１であるならば、Ｐ１Ｌは単に単位行列（identity matrix）とすることができることが理解されよう。これは、フレームＲ１が、１つの軸に沿って１つ（または複数）の自由度を表現するのに適する場合であろう。
【０１２４】
したがって、ハイブリッドリファレンスフレームを用いて、式（３）および（４）は、それぞれ、
【０１２５】
【数１５】
2006294038.tif 000017
、および、
【０１２６】
【数１６】
2006294038.tif 000018
【０１２７】
と表される。これにより、２つの有限要素に関する以下の剛性行列が導き出される。
【０１２８】
【数１７】
2006294038.tif 000019
【０１２９】
式（８）は、制約が緩和されていない場合において用いられる計算と実質的に同一の計算により、剛性行列と、システムに与えた力とが、ハイブリッドリファレンスフレームを用いて表現できることを示している。わずかな違いは、ローカルリファレンスフレームＲＬを見つけなければならないこと、および、ローカルリファレンスフレームＲ１およびＲ２を、制約が緩和されている、部分的に制約された節点（群）における新たなローカルリファレンスフレームに変換するためのさらなる２つの行列Ｐ１ＬおよびＰ２Ｌを計算しなければならないことである。
【０１３０】
しかしながら、式（８）が計算されると、エネルギを最小化するのに、上記にて与えられた式（２）における行列ＮおよびＭの計算と同様の計算は必要ない。実際、緩和された自由度がリファレンスフレームの１つの軸に沿って表現されるように、ローカルリファレンスフレームＲＬは選択される。これは、サブセットＸａの選択が、緩和されているすべての制約を単に選択することによって、実行され得ることを意味する。すなわち、式（８）において、制約が緩和されていないローカルリファレンスフレームの１つまたは複数の軸に対応する、剛性行列の行および列と、荷重ベクトルの行とが単に取り除かれる。
【０１３１】
具体的には、図６の例において、各節点ごとに、３つの並進自由度（ｕ，ｖ，ｗ）と、３つの回転自由度（θｘ，θｙ，θｚ）とからなる６自由度を考慮する。各Ｘｉは、以下のように記述される。
【０１３２】
【数１８】
2006294038.tif 000020
【０１３３】
したがって、ハイブリッドリファレンスフレームを用いて表現された、式（８）における剛性行列は、２４×２４の行列である。節点４１において１つの自由度が存在し、節点は、軸ｕのまわりを自由に回転することができる、すなわち、θｘ２は緩和されていると仮定されたい。これは、ｕ２、ｖ２、ｗ２、θｙ２およびθｚ２は固定されているのに対し、θｘ２と、Ｘ１およびＸ３におけるすべての変数とは計算される必要があることを意味する。端節点を除いた節点における変位を計算するため、このような趣旨で、図２および図３を参照して上記にて開示したように、処理を進める。
【０１３４】
より具体的には、まず、式（８）が、以下のように詳細に記述される。
【０１３５】
【数１９】
2006294038.tif 000021
【０１３６】
この式において、１≦ｉ≦１８かつ１≦ｊ≦１８であるｋｉｊは、上述したように計算された、剛性行列の様々な係数である。ｕ、ｖ、ｗ、θｘ、θｙおよびθｚに対する１番目の添え字は、節点を表している。ｕ、ｖ、ｗ、θｘ、θｙおよびθｚに対する２番目の添え字は、変位が表現されるフレームを表していて、Ｇは、グローバルリファレンスフレームＲＧを表し、Ｌは、ローカルリファレンスフレームＲＬを表す。
【０１３７】
印加荷重に関しては、Ｆｘａｐｐ、ＦｙａｐｐおよびＦｚａｐｐが、節点に与えた力のコンポーネントを表しているのに対し、Ｍｘａｐｐ、ＭｙａｐｐおよびＭｚａｐｐは、節点に与えたトルクのコンポーネントを表している。１番目の添え字（アラビア語の添え字）は、荷重が与えられた節点を表している。２番目の添え字は、変位が表現されるフレームを表していて、Ｇは、グローバルリファレンスフレームＲＧを表し、Ｌは、ローカルリファレンスフレームＲＬを表す。
【０１３８】
局所座標系ＲＬを用いて表現された５つの制約条件
【０１３９】
【数２０】
2006294038.tif 000022
【０１４０】
と併せて考慮されると、式は以下のようになる。
【０１４１】
【数２１】
2006294038.tif 000023
【０１４２】
エネルギが最小化される（したがって、中央の節点におけるローカルなｕ軸についての緩和された回転、すなわち、ローカルな自由度θｘ２を、他の２つの節点におけるすべてのグローバルな自由度とともに含む、自由度のＸａセットを得る）と、以下のように自由度のフルベクトル（full vector）Ｘを復元することができる。
【０１４３】
【数２２】
2006294038.tif 000024
【０１４４】
この復元式において、変位ベクトルコンポーネントが、部分的にはローカルリファレンスフレームを用いて表現され、すなわち、中央の節点において表現され、部分的にはグローバルリファレンスフレームを用いて表現されている、すなわち、他の２つの節点において表現されていることに留意されたい。ここで、グローバルリファレンスフレームを用いて完全に変位ベクトルを表現することが可能であることは言うまでもない。
【０１４５】
Ｘａベクトルのすべてのコンポーネントが計算されると、必要に応じて、上述したように、共回転アップデートに加えて節点リファレンスフレームのアップデートを実行することができる。唯一の違いは、自由度を有する節点に対して、アップデートが、計算されたコンポーネントと、制約されたコンポーネントとの組合せを用いることである。図６の例では、ｕ２、ｖ２、ｗ２、θｙ２およびθｚ２が固定されているのに対し、θｘ２は、Ｘ１およびＸ２のコンポーネントと同時に計算される。リファレンスフレームのアップデートおよび節点フレームのアップデートに関しては、変位コンポーネント（ｕ２，ｖ２，ｗ２，θｘ２，θｙ２，θｚ２）が、
−予め定められた変位境界条件の一部と、
−インクリメンタルステップ中のエネルギの最小化と
から生じるという事実を考慮することなく、これらのコンポーネントを用いて、上記にて開示したように、処理を進めることができる。アップデートが必要とされない場合、例えば、線形システムの場合または連続体システムの場合、次のステップに直接進むことができる。
【０１４６】
本発明によって、システムの有限要素表現の制約された節点におけるいくつかの自由度を緩和することが可能となる。これは、エネルギが最小化されるフレームを慎重に選択することによって、制約のあらゆるファクトリゼーションを回避することが可能となる、という事実を用いる。
【０１４７】
図７ないし図２２を参照しながら、以下にプロセスの実施例を提供する。これら４つの実施例は、フランスのＳｕｒｅｓｎｅｓを本拠とするＤａｓｓａｕｌｔ  Ｓｙｓｔｅｍｅｓ社のＣＡＴＩＡ（登録商標）  Ｖ５内に統合された非線形ソルバコア（non-linear solver core）（ＦＯＲＴＲＡＮコード）を用いる。各実施例について、シミュレーションは、７５回のインクリメンタルステップ（開始コンフィギュレーションから初期コンフィギュレーションまで４５回のステップ、初期コンフィギュレーションから最終コンフィギュレーションまで３０回のステップ）からなる。これらのステップのうちの少数のステップのみが図に示されている。すべての処理は、ＩＢＭ（登録商標）  Ｔｈｉｎｋｐａｄ（登録商標）  Ｔ４２（１．８ＧＨｚ  Ｐｅｎｔｉｕｍ  ＩＶ（登録商標）プロセッサ、１．０ＧＢのＲＡＭ）上で実行された。
【０１４８】
第１の実施例は、図７ないし図１０に示されている。システムは１つのセグメントである。両端点は完全に制約される。システムは１２の要素モデルとして解析される。完全に処理を実行するのに０．２８秒のＣＰＵ時間を要する。初期段階から最終段階に至る各インクリメンタルステップにおいて、平均６回の反復をもって収束が達成される。
【０１４９】
節点の位置および節点リファレンスフレームが、すべての図において示されている。図７は開始コンフィギュレーションを示しており、この場合、システムの形状は、ユーザによって完全に指定される。
【０１５０】
図８は初期コンフィギュレーションを示している。端点の位置および向きは、ユーザによって指定され、図では、薄いリファレンスフレームとして示されている。システムの変形した形状は、プログラムによって計算される。
【０１５１】
図９は最終コンフィギュレーションを示している。端点の位置および向きは、ユーザによって指定され、濃いリファレンスフレームとして示されている。変形した形状は、プログラムによって計算される。図８における端点の位置および向きと、図９における端点の位置および向きとの違いは、システムが受ける予め定められた変位境界条件を表す。
【０１５２】
図１０は３つの中間コンフィギュレーションを示している。端点の位置および向きが補間され、この位置および向きは、予め定められた変位境界条件の一部を表す。変形した形状は、プログラムによって計算される。
【０１５３】
第２の実施例は、図１１ないし図１４に示されている。初期コンフィギュレーションおよび最終コンフィギュレーションの双方において、（変形平面に直交するローカルな軸まわりの）回転自由度が左端点にて解放されていることを除き、これは第１の実施例と類似している。
【０１５４】
完全に処理を実行するのに０．２７秒のＣＰＵ時間を要する。第１の実施例と同様に、平均して、インクリメントにつき６回の反復をもって収束が達成される。
【０１５５】
図１１−図１４は、それぞれ、図７−図１０に対応するものであり、これらについては再度説明しない。第１の実施例に関しては、節点の位置および節点リファレンスフレームが、すべての図において示されている。
【０１５６】
第３の実施例は、図１５ないし図１８に示されている。この実施例では、システムは、異なる特性の３つのセグメントから構成される。１つの端点、および、１つの中間点が、完全に制約される。これは図に表され、この図において、システムのより細長い部分は、制約されておらず、重力の作用下で、下方に曲がっていることがはっきりと分かる。
【０１５７】
有限要素解析が３×１２の要素モデルに対して実行される。
【０１５８】
完全に処理を実行するのに６．４５秒のＣＰＵ時間を要する。平均して、インクリメントにつき４０回の反復をもって収束が達成される。
【０１５９】
図１５−図１８は、それぞれ、図７−図１０に対応するものであり、これらについては再度説明しない。節点の位置および節点リファレンスフレームが、すべての図において示されているわけではない。
【０１６０】
第４の実施例は、図１９ないし図２２に示されている。３回転自由度のすべてが中間点にて解放されていることを除き、この実施例は第３の実施例と同一である。
【０１６１】
完全に処理を実行するのに５．２秒のＣＰＵ時間を要する。各インクリメンタルステップにおいて、平均して、インクリメントにつき２５回の反復をもって収束が達成される。
【０１６２】
図１９−図２２は、それぞれ、図１５−図１８に対応するものであり、これらについては再度説明しない。システムに与えた異なる制約に起因する、最終コンフィギュレーションにおける違いに留意されたい。
【０１６３】
プロセスは、上記にて提供した実施例に限定されるものではない。例えば、図１の例では、予め定められた変位が、フレキシブルシステムの両端において考慮される。システム内の３以上の点において、予め定められた変位（完全に制約されている、または、制約されていない）を考慮してもよい。
【０１６４】
図１ないし図３を参照して開示したプロセスにおいて、共回転アップデートは、各インクリメンタルステップの後に、すなわち、次のインクリメンタルステップの前に、実行される。各ステップにおいて共回転アップデートを実行しないことも可能である。具体的には、必要なときだけ、すなわち、各エレメンタリリファレンスフレームと、対応する最新のエレメンタリリファレンスフレームとの間の角度の変化が予め定められた値よりも大きいときだけ、共回転アップデートを実行してもよい。事実上、この予め定められた値は、小さな値、例えば１秒に対して、０度以上３０度以下とすることができ、これは、擬似的なシステマティックアップデーティング（quasi-systematic updating）に対応する。
【０１６５】
インクリメンタルステップ中の収束を制御するためのスケーリング係数が、上記にて開示されている。このようなスケーリング係数は、フレキシブルシステムの挙動を制御するために使用することもできる。例えば、所与の一インクリメンタルステップにおいてゼロの値をスケーリング係数に適用すると、粘弾性フレキシブルシステムをシミュレートする結果となる。すべての初期応力が緩和されている。
【０１６６】
フレキシブルシステムのより効率的なシミュレーションのために、上記にて開示した様々な実施形態を組み合わせることができる。様々な実施形態を個別に実行することもできる。例えば、収束判定条件は、節点フレームとは無関係に使用することができる。上述した収束判定条件は、図３および図４を参照して説明した共回転アップデートとは無関係に使用することができる。収束判定条件はまた、非構造要素（non-structural element）（すなわち、連続要素）に対して使用することもできる。このような非構造要素、すなわち連続要素は、実際には「ソリッド」要素であり、節点の回転自由度を有さない。
【０１６７】
上述したように、図６を参照して開示したハイブリッドリファレンスフレームは、前述した共回転アップデートとは無関係に使用することができる。このようなハイブリッドリファレンスフレームはまた、線形ソルバにおいて使用することもできる。
【図面の簡単な説明】
【０１６８】
【図１】開始コンフィギュレーションおよび最終コンフィギュレーションにおけるフレキシブルシステムの概略図であって、複数の中間コンフィギュレーションを示す図である。
【図２】開始コンフィギュレーションおよび第１中間コンフィギュレーションにおける、フレキシブルシステムに合わせた有限要素モデルの部分図表現である。
【図３】図２の有限要素の拡大図である。
【図４】様々なリファレンスフレームを伴う有限要素の図である。
【図５】様々なリファレンスフレームを伴う有限要素の図である。
【図６】システムの節点が自由度を有する場合に使用されるフレームを伴う、２つの有限要素の図である。
【図７】プロセスの実施例を示す図である。
【図８】プロセスの実施例を示す図である。
【図９】プロセスの実施例を示す図である。
【図１０】プロセスの実施例を示す図である。
【図１１】プロセスの実施例を示す図である。
【図１２】プロセスの実施例を示す図である。
【図１３】プロセスの実施例を示す図である。
【図１４】プロセスの実施例を示す図である。
【図１５】プロセスの実施例を示す図である。
【図１６】プロセスの実施例を示す図である。
【図１７】プロセスの実施例を示す図である。
【図１８】プロセスの実施例を示す図である。
【図１９】プロセスの実施例を示す図である。
【図２０】プロセスの実施例を示す図である。
【図２１】プロセスの実施例を示す図である。
【図２２】プロセスの実施例を示す図である。
【図１】
2006294038.tif 000025
【図２】
2006294038.tif 000026
【図３】
2006294038.tif 000027
【図４】
2006294038.tif 000028
【図５】
2006294038.tif 000029
【図６】
2006294038.tif 000030
【図７】
2006294038.tif 000031
【図８】
2006294038.tif 000032
【図９】
2006294038.tif 000033
【図１０】
2006294038.tif 000034
【図１１】
2006294038.tif 000035
【図１２】
2006294038.tif 000036
【図１３】
2006294038.tif 000037
【図１４】
2006294038.tif 000038
【図１５】
2006294038.tif 000039
【図１６】
2006294038.tif 000040
【図１７】
2006294038.tif 000041
【図１８】
2006294038.tif 000042
【図１９】
2006294038.tif 000043
【図２０】
2006294038.tif 000044
【図２１】
2006294038.tif 000045
【図２２】
2006294038.tif 000046
【外国語明細書】
2006294038000001.pdf