JP4653482B2 - コンピュータ変形解析器 - Google Patents

Description

この発明は概してコンピュータ支援によるモデリングに関わるものであり、更に詳しくは、コンピュータ変形解析器（ｃｏｍｐｕｔｅｒｉｚｅｄ ｄｅｆｏｒｍａｔｉｏｎ ａｎａｌｙｚｅｒ）に関わるものである。

機械装置は多くの場合、有限要素解析（「ＦＥＡ」）パッケージを用いて設計される。例えば、金属薄板を覆うメッシュを生成することで、このＦＥＡパッケージは自動車のフード（ｆｏｏｄ）のモデル形成に使用される。その後ＦＥＡパッケージは、メッシュ全体の結果として得られる形がはっきりと描かれるようになるまで、ある境界条件に基づきメッシュの各要素の偏移を徐々に決定する。モデルが導く結果は実際の生産プロセスで使用されるので、ＦＥＡパッケージによって生成されるこのモデルの精度は重要である。例えば、モデルが正確でなければ、そのモデルに基づいて生産された実際のデバイスは、好ましい形状を成さない可能性がある。

ＦＥＡパッケージで生成されたモデルの精度は、モデリングプロセスで使用された特定のメッシュを形成する要素の数を増やすことで向上する。しかし、要素数の増加により、演算処理所要時間も増加する。何故なら、非直線形のＦＥＡではメッシュの各要素に関連した多くの式の演算を繰り返し行うからである。これらのモデルの演算時間を短縮するためにスーパーコンピュータが使用されるが、こういった強力なコンピュータでさえ、要求精度を有するモデルを完成させるには数日かかる。

本発明の一つの実施例には、素材表面の変形を決定する方法が示されている。この方法には、コンピュータによる素材の弾性変形範囲及び塑性変形範囲の受信が含まれる。この方法には、メッシュでもって素材の領域をモデリングすることも含まれる。メッシュには複数の要素が含まれる。この方法にはまた、コンピュータを使用しメッシュの有限要素解析を行うことで、各要素の少なくとも一つのポイントの変位を決定することも含まれる。この変位はモデリングされた一連の境界条件から発生する。有限要素解析は、弾性変形範囲及び塑性変形範囲に基づき各要素の弾性係数を反復的に近似することを含む。近似化は、事前に計算された要素の変位から計算された歪みの値に基づき求められる。

本発明の実施例の幾つかは、多数の技術的効果を提供する。実施例の中には、これらの効果の幾つか、または全ての恩恵を受けているものもあり、その恩恵を受けていないものもある。例えば、幾つかの実施例によると変位の大幅な増加があるモデルでさえ、メッシュの要素数を大幅に増加させること無く、デバイスの正確なモデルを生成できる。あるメッシュの要素数が比較的小さければ、要求される演算の能力と演算時間が少なくてすむ。別の実施例によると、より高位の要素を使用することによっても、多数の要素の必要性が減るので、正確なモデリングに必要な演算の能力と演算時間が少なくてすむ。
その他の技術的効果は、当業者には簡単に確認できる。

本発明の実施例は、図１Ａ乃至図５Ｄの図面を参照することでよく理解できる。それぞれの図で同一で対応する部分には同一の参照番号を使っている。

図１Ａ及び図１Ｂは、図示する素材が特定の形状への形成をどうモデリングするかを図示する。素材の一例として、金属薄板（ｓｈｅｅｔ ｍｅｔａｌ）がある。図１Ａに示すメッシュ（ｍｅｓｈ）５０Ａは、形成過程前の素材のモデルの実施例である。図１Ｂに示すメッシュ５０Ｂは、形成過程後の素材のモデルの実施例である。

メッシュ５０Ａ及び５０Ｂは、有限要素解析（ＦＥＡ）のソフトウェアプログラムまたは関連するメッシュジェネレータ（明示せず）で生成できる。メッシュ５０Ａ及び５０Ｂは、パラソリッドファセット形成を含む如何なる適切なプロセスをもってしても形成できる。図１Ａに示す実施例においては、メッシュ５０Ａは４つの角５２を有するが、メッシュのジオメトリー次第で、角５２の数はこれ以上にも以下にもなりえる。

メッシュ５０Ａは、これから変形する素材の特定の領域５４のモデルである。領域５４は、目的の形状を有するよう変位する素材の部分である。メッシュ５０Ａが素材の領域５４の変位のモデルとなるようにメッシュ５０Ａの一部分を適切に変位させるために、一つ以上の境界条件がメッシュ５０Ａに課すことができる。メッシュ５０Ａは複数の要素５８から成る。要素５８は、三角形、四角形、またはその他の多角形などの簡単なジオメトリーを有することができる。図１Ａのメッシュ５０Ａには四角形のみが示されているが、素材の特定のジオメトリーのモデルとなるメッシュ５０Ａを形成するために、如何なる多角形、または多角形の組み合わせも使用することができる。各要素５８は一つ以上のガウスポイント６４から成る。このガウスポイントは、符号のインテグレーションのために選択されたポジションのことである。特定の場合のガウスポイント６４の数は、要求されるモデルの精度及び要素５８のジオメトリーに基づき当業者なら選択できる。各要素５８は少なくとも一つの隣接した要素の端部６０により他の要素５８に繋がっている。

図１Ｂを参照すると、メッシュ５０Ｂは、初期のメッシュ５０Ａの４箇所の角５２に固定させ、中心５３を上方向に変位させて形成した結果得られるものである。角５２の固定及び中心５３の変位はＦＥＡプロセスの境界条件として課すことができる。上方向に変位するのは中心５３のみであるが、結果としてこの変位は中心５３の周りの素材表面及び素材の端部を変位させる。これらの変位は、角５２を固定しているにも拘わらず中心５３を動かしたことに対する反応である。素材が異なれば、発生する変位のタイプもレベルも異なる。例えば、ゴムでは端部及び動かされた特定の領域を取り囲む部分において大きな変位が発生し得る。何故なら、ゴムは永久塑性変形素材の特徴を持った挙動ではなく、ゴム素材の特徴を持った挙動を取るからである。また、メッシュ５０Ａに異なる制限と動きが課されるよう、異なる境界条件を課すこともできる。異なる一連の境界条件を与えることで、結果として得られる最終的なメッシュの形状が図２Ｂに示されるメッシュ５０Ｂのそれとは異なることになる。

曲げ、伸展といった異なる素材の異なる動きの最終的な変位を正確に予測するために、従来のＦＥＡ法では何千もの要素を有するメッシュを使用する。本発明の内容から分かることは、このような多大な要素数は、ＦＥＡプロセスによって複製されていない条件によって生じ得る最終的なモデルの不正確さを補うために必要だということである。複製されていない条件の例としては、素材特性の変化、特に徐々に変位が大きくなることによる弾性係数の変化が挙げられる。他の従来型ＦＥＡパッケージの多くは小規模の変位にのみ有用であり、そのため、素材の特性を変化させるような大きな変位を分析するツールとしては精度に欠ける。更にまた、多くの要素を有するということは、演算時間と演算能力の大幅な増加を必要とする。このような増加は避けられないものである。何故なら、各要素は少なくとも一つのシェープ関数と言われる、素材の挙動を数学的にモデリングするために使用される数式と関連しているからである。これら全ての式は、全式が一つの解に収束するまで反復的に一緒に解かれる。製造業界では要求される演算能力を満たすために、十分な精度を有するデバイスのモデルが入手できるまで繰り返しこれらの演算を行うために、しばしばスーパーコンピュータを使用する。しかしながら、スーパーコンピュータを使用した場合でも、適切なモデルを入手するには数日間の演算時間が必要となる。

本発明の実施例の幾つかでは、イテレーション的なＦＥＡプロセス中に各要素内で弾性係数などの素材の不均一な特性をアップデートすることで、デバイスのモデルの精度を向上させる方法を示している。大きな変位が発生するモデルにおいてでさえ、メッシュ内の要素数を大幅に増加させずにデバイスの正確なモデルが作成できると言う点で、このことには優位性がある。あるメッシュ内で要素数が比較的少ないと、必要な演算能力と演算時間は減少する。システムの実施例及び方法の実施例の更なる詳細を、図２乃至図５Ｄと合わせて以下に更に詳細に説明する。

図２は、例としての素材の応力と歪みとの関係の観点から、弾性係数を図示した例としてのグラフ９０である。「弾性係数」という言葉は、弾性変形範囲及び塑性変形範囲両者内での、素材の応力と歪みの比率を意味する。「弾性変形範囲」または「弾性ドメイン」とは、固形物の形状や大きさの変化が可逆的である歪みの範囲を意味する。「塑性変形範囲」または「塑性ドメイン」とは、弾性ドメインの応力限度を超える歪み範囲を意味する。「塑性ドメイン」つまり「塑性変形範囲」内での固形物の形状または大きさの変化は恒久的なものであるが、その固形物の破損は発生しない。軸９４は応力（「σ」）を示し、軸９８は、歪み（「ε」）を示す。形成プロセス中に大きな変位が発生し得る実施例の幾つかにおいては、素材の弾性ドメインから塑性ドメインへの移行を説明する必要がある場合もある。実施例の一つにおいては、素材の降伏応力（「σ_y」）を示すポイント１０８は、素材の弾性ドメインを示すセグメント１００から素材の塑性ドメインを示すセグメント１０４への移行を示している。一つの実施例においては、この弾性―塑性関係の近似化に、弾性ドメインを線形の応力―歪み関係で近似し、塑性ドメインを別の線形な応力―歪み関係で近似する、準線形アプローチが適用される。弾性ドメインであるセグメント１００の傾きは、弾性係数（「Ｅ」）で表される。塑性ドメインセグメント１０４の傾きは、接線係数（「Ｅ_Ｔ」）で表される。

グラフ９０に示されるモデルによると、応力―歪み関係は以下のように定義される。

ここでε_ＰＬは、ポイント１０８で歪みと等しくなる。

と決定できる。ここでＨ’＝有効応力―歪みの傾き、である。

図２に示され、また式１、２に定義される応力―歪み関係などの、素材特性を提示することは有益である。何故ならその結果得られるモデルが、素材の実際の変形した形状をより正確に反映するからである。それに加え、応力―歪み関係を準線形関係として近似させることは有益である。何故なら、弾性変形と塑性変形のこのような線形近似化はＦＥＡに含まれる計算を簡素化するからである。しかしながら、素材の応力―歪み関係は非線形関数によって示されることも、実験データによって示されることもあり得る。数学的簡略化により、要求される演算能力と演算時間が少なくてすむ。

図３Ａは、本発明の内容の恩恵を受けた有限要素解析を行うためのコンピュータシステム１０の実施例を図示した概略図である。システム１０は一つ以上の入力デバイス１８及び、一つ以上の出力デバイス２０を接続したコンピュータ１４を含む。ユーザー２４はシステム１０にアクセスでき、データを入力するためや、図面２８を作成編集するために入力デバイス１８を使用でき、図面２８は出力デバイス２０の何れか、または全てに表示される。

図３Ａに示されるように、入力デバイス１８の例として、キーボードとマウスが挙げられる。しかしながら、入力デバイス１８は、スタイラス、スキャナ、またはそれらの何れかの組み合わせを含め、他の形態でもよい。出力デバイス２０の例としては、色々なタイプのモニターとプリンタが挙げられる。しかしながら、出力デバイス１６は、プロッタを含め他の形態をとることもできる。液晶ディスプレイ（「ＬＣＤ」）、またはブラウン管ディスプレイ（「ＣＲＴ」）など、ユーザー２４が図面２８を見ることができる適切なビジュアルディスプレイユニットはどんなものでも出力デバイス２０として適切である。

図３Ｂは本発明の一つの実施例による有限要素解析を行う上で使用されるコンピュータ１４のブロック構成図である。図示されているように、コンピュータ１４はプロセッサ３０、ＦＥＡソフトウェアプログラム３８を記憶するメモリ３４、及びＦＥＡソフトウェアプログラム３８に関連するデータまたは他のデータ記憶する一つ以上のデータ記憶ユニット４０を含む。

プロセッサ３０は、メモリ３４及びデータ記憶ユニット４０に接続されている。プロセッサ３０は、ＦＥＡソフトウェアプログラム３８のロジックを実行するために、また図面に関するデータを呼び出すべくまたは記憶すべくデータ記憶ユニットにアクセスするために作動可能である。プロセッサ３０の例として、インテル社が市販しているペンティアム（登録商標）^ＴＭシリーズのプロセッサが挙げられる。

メモリ３４及びデータ記憶ユニット４０はファイル、スタック、データベースまたは他の適切なデータ形式から成ってもよい。メモリ３４、及びデータ記憶ユニット４０はランダムアクセスメモリ、読み取り専用メモリ、ＣＤ−ＲＯＭ、取り外し可能なメモリデバイス、あるいはデータの記憶及び／または取り出しが可能な適当なデバイスであればどのようなものでもよい。メモリ３４とデータ記憶ユニット４０は入れ替えてもよく、同一の機能を果たすことができる。

ＦＥＡソフトウェアプログラム３８は、ユーザー１８がコンピュータ１４を使用して有限要素解析でデバイスのモデリングを行うことを可能にするコンピュータプログラムである。ＦＥＡソフトウェアプログラム３８は、コンピュータ支援製図（「ＣＡＤ」）パッケージなどの製図アプリケーションの一部であってもよいし、または独立したコンピュータアプリケーションであってもよい。ＦＥＡソフトウェアプログラム３８は、例えばメモリ３４やデータ記憶ユニット４０のような、如何なる記憶媒体内に存在してもよい。ＦＥＡソフトウェアプログラム３８は、適切なコンピュータ言語であればＣやＣ＋＋を含むどのようなコンピュータ言語で書かれていてもよい。ＦＥＡソフトウェアプログラム３８は、ユーザー１８が曲線の初期位置や曲線の最終位置などの境界条件を入力し、その結果得られる予想形状を出力デバイス２０に表示し及び／またはデータ記憶ユニット４０に記憶するためにするために作動可能である。本発明の内容を取り入れることのできるＦＥＡソフトウェアプログラム３８の例を挙げると、ユニグラフィックスソリューションズ社が市販しているリージョンアナライザー^ＴＭがある。

図４Ａは、図２に示される初期のメッシュ５０Ａのようなメッシュの四角形の要素１３０Ａの実施例の概略図である。幾つかの実施例においては、要素１３０Ａはシェル要素でもよい。シェル要素とは、図４に示されるように、それと関連する厚さを有するジオメトリーを意味する。要素１３０Ａは、要素１３０Ａの要素の端部１６０、１６４、１６８、及び１７０に位置する複数のノード１３２乃至１５４から成る。ノード１３２乃至１５４は、まとめてノード１３１と呼ぶ。ノード１３１は、要素１３０Ａの様々な部位の変位を決定する際の基準として使用される。図４Ａに図示する実施例においては、要素１３０Ａは１２のノード１３１を有するが、これより多い数のノード、または少ない数のノードを使用してもよい。ノード１３２、１３８、１４４、及び１５０は角のノードである。ノード１３４、１３６、１４０、１４２、１４６、１４８、１５２、及び１５４は辺の途中のノードである。ノード１３１のローカル座標の位置を、要素１３０Ａのモデル１３０Ｂに示す。

ＦＥＡプロセスの一部として、ノード１３１のそれぞれのシェープ関数（ｓｈａｐｅ ｆｕｎｃｔｉｏｎ）Ｎが決められる。ノード１３１に関連するシェープ関数を以下に説明する。これらの式で使われているように、図４Ａにおける参照番号とノード番号（i）の間には、１３２＝ノード１、１３４＝ノード２、１３６＝ノード３、１３８＝ノード４、１４０＝ノード５、１４２＝ノード６、１４４＝ノード７、１４６＝ノード８、１４８＝ノード９、１５０＝ノード１０、１５２＝ノード１１、及び１５４＝ノード１２、という関係がある。

図４Ｂは、図２に示される初期のメッシュ５０Ａのような初期のメッシュの三角形の要素１９０Ａの実施例の概略図である。幾つかの実施例においては、図４Ｂに示されるように要素１９０Ａはシェル要素である。要素１９０Ａは、要素１９０Ａのノードの線２１２、２１４、及び２１８に位置する複数のノード１９２乃至２０８から成る。ノード１９２乃至２０８は、まとめてノード１９１と呼ぶ。ノード１９１は、要素１３０Ａの様々な部位の変位を決定する際の基準として使用される。図４Ｂに示される実施例においては、要素１９０Ａは１０のノード１９１を有するが、これより多い数のノード、または少ない数のノードを使用してもよい。ノード１９２、１９８、及び２０４は角のノードである。ノード１９４、１９６、２００、２０２、２０６、及び２０８は辺の途中のノードである。ノード２１０は、中心のノードである。ノード１９２乃至２１０のローカル座標を、要素１９０Ａのモデル１９０Ｂに示す。ＦＥＡプロセスの一部としてノード１９１のそれぞれのシェープ関数Ｎが決定され、それらを以下に記述する。これらの式に使われているように、参照番号とノード番号（i）の間には、１９２＝ノード１、１９４＝ノード２、１９６＝ノード３、１９８＝ノード４、２００＝ノード５、２０２＝ノード６、２０４＝ノード７、２０６＝ノード８、２０８＝ノード９、及び２１０＝ノード１０、という関係がある。

ここで、Ｎ_ｉ、Ｌ_１、Ｌ_２、Ｌ_３は、次の通りである。
Ｎ_ｉ＝ノードｉに関連するシェープ関数。
Ｌ_１＝第１ノードと隣接する端部間のシェープ関数の大きさ。
Ｌ_２＝第２ノードと隣接する端部間のシェープ関数の大きさ。
Ｌ_３＝第３ノードと隣接する端部間のシェープ関数の大きさ。

幾つかの実施例においては、１０のノードを有する三角形の要素１９０Ａのシェープ関数を、３つの独立変数（Ｌ_１、Ｌ_２、Ｌ_３）ではなく、２次元のパラメータ領域を表すし２つの独立変数（ｕ、ｖ）で定義することが望ましい。Ｌ_１＋Ｌ_２＋Ｌ_３＝１という関係を使い、これらのシェープ関数を独立変数Ｌ_２、Ｌ_３に関して次のように書き換えることもできる。

ここで、Ｌ_１を１−Ｌ_２−Ｌ_３に置き換えた。モデル１９０Ｃは（ｕ、ｖ）のフォーマットでのノード１９２乃至２１０を示す。
幾つかの実施例においては、Ｕが要素のＵ方向のパラメータの値を示し、Ｖが要素のＶ方向のパラメータの値を示すとして、Ｕ＝Ｌ_２、Ｖ＝Ｌ_３という観点からこれらの式を、次の通りのように書き換えると便利であることもある。

一次偏導関数は以下のように計算できる。

要素１３０Ａは要素の例として使用されているが、他の形状、ノード数、または厚さを有する要素もまた使用可能である。例えば、要素１９０Ａを、要素１３０Ａと組み合わせて使ってもよい。

図５Ａは本発明の一つの実施例によるＦＥＡの方法２５０を図示するフローチャートである。この方法２５０は、有限要素ソフトウェアプログラム３８または他の適切なソフトウェアまたはハードウェアで実行できる。方法２５０はステップ２５４から開始する。ステップ２５８では、以下グローバル変位マトリクス｛Ｄ｝と呼ぶグローバル変位を使用し、ローカル要素座標が確立される。ローカル要素座標は、要素のローカルＵＶ方向と、要素の境界のどのポジションであろうと、そのポジションでの法線ベクトルに沿った座標システムのことを意味する。グローバル変位とは、単一の「グローバル」座標システムに対する全ノードの変位のことを意味する。ステップ２５４の更なる詳細を、図５Ｂと合わせ後に示す。

ベクトル｛ｄ｝で表す要素の変位は、適切な数学的プロセスを使いステップ２６０でローカル座標で計算される。要素の変位の計算は、ステップ２５８で計算されるグローバル変位と、ステップ２５８で計算され確立されたローカル座標に基づいている。

ステップ２６４では、剛性マトリクス［ｋ］と、計算された要素の変位｛ｆ｝＝−［ｋ］｛ｄ｝から結果として生じる力の両者が、ローカル座標で計算される。ステップ２６４の更なる詳細を図５Ｃと合わせ後に記述する。これらの変換されたバージョンは、不正確に変位を近似した結果発生する要素５８の不均衡力に基づき、ノード１３１の変位の次のイテレーションを決定するためにその後使用される。

ステップ２６８では、それぞれの要素５８の要素剛性マトリクス［ｋ］と力ベクトル｛ｆ｝がグローバル座標に変換される。ステップ２７０では、プロセスされていない要素１３０Ａがあるかどうかを確認する。プロセスされていない要素１３０Ａが残っている場合、方法２５０のステップ２５８へ戻る。プロセスされていないそれぞれの要素１３０Ａに関して、ステップ２５８乃至２６８がステップ２７０でプロセスされていない要素が無いと確認されるまで繰り返される。

ステップ２８４において、変位の増加分｛ΔＤ｝はグローバル変位｛Ｄ｝に加えられる。ステップ２８６では、収束が達成されたかを決定するために、ノード１３１各々の結果として得られた変位の増加分を、それぞれの対応するノード１３１の直前の変位と比較する。一つの実施例においては、全てのノード１３１のその時点での変位と以前の変位との差が事前に決められたモデリングの許容範囲よりも小さければ、収束が達成し得る。モデリングの許容範囲は、一般的にゼロに近い値に設定される。例えば、一つの実施例ではモデリングの許容範囲は０．００１インチに設定されるが、適切な許容範囲は如何なるものであれ適用できる。収束が達成されていない場合は、ステップ２５８乃至２８６が繰り返される。このようにして、下記のイテレーション的な関係が確立する。

ここで、ｉはイテレーションの回数を表す。

ステップ２５８乃至２８６のそれぞれのイテレーションにおいて、各要素１３０Ａのそれぞれの剛性マトリクス［ｋ］は、ステップ２６４で再度評価される。要素１３０Ａの剛性マトリクス［ｋ］は、ステップ２５８乃至２８６のそれぞれのイテレーション後の要素５８の変位レベルによって変化し得る。変位の増加分と共に、要素５８の各ガウスポイントでの要素の剛性は、要素のノード変位によって発生する実際の歪みの量により変化し得る。例えば、降伏応力を超してしまい、そのため素材の剛性が弾性素材の挙動から塑性素材の挙動へと移行することもあり得る。ステップ２５８乃至２８６のイテレーション後、各ガウスポイントに、そのガウスポイントが弾性ドメインと塑性ドメインのどちらにあるかに基づき、新しい剛性関係を割り当てることにより、要素５８の各ガウスポイントでの剛性影響を計算し、各要素の最新の剛性マトリクスを決定することができる。この代わりとして、素材剛性の新しい数値をイテレーション毎よりも少ない頻度でアップデートすることもできる。ステップ２５８乃至２８６のイテレーション毎に、後で記述する図５Ｂのグラフ７０が示すような素材の応力―歪み関係を用いて、弾性係数をアップデートできる。ステップ２５８乃至２８６のイテレーションにおける弾性係数及び剛性マトリクスのアップデートは、弾性係数及び剛性マトリクスの「反復的近似化」または「反復的アップデート」と呼ぶこともできる。

ステップ２５８乃至２８６のイテレーションのために剛性マトリクス［ｋ］をアップデートすることは有益である。何故ならこれが、特定の素材の弾性ドメインと塑性ドメインとの間の移行を補完するからであり、この例を図２に示す。このような補完は、変位の増加分が大きい場合でさえ、メッシュ中の要素数を増加させることなく更に精度の高いモデルを結果として生み出す。

ステップ２８６で収束が達成されたと決定されれば、ステップ２８８で弾性のスプリングバック／オーバークラウン解析をオプションとして行うことができる。ステップ２８８の更なる詳細を後に説明する。方法２５０はステップ２４０で完了する。

ステップ２６４のイテレーション時に剛性マトリクス［ｋ］を再評価及びアップデートすることは有益である。何故なら、あるメッシュの要素数を大幅に増加させることなく、正確なデバイスのモデルが生成できるからである。あるメッシュ中の要素数が比較的少ない場合は、要求される演算能力と演算時間が少なくてすむ。

図５Ｂは、要素５８のような要素の弾性係数の放物線近似を図示するグラフ７０である。一つの実施例においては、要素には各ガウスポイント６４に割り当てられた弾性係数の独立な値がある（図１Ａ及び図１Ｂに示す）。要素の変位前では、各ガウスポイント６４での弾性係数の大きさの分布は、平らになっているグラフ７０の点線部分に示されるように、平坦で平らなものとして図示されている。増加分だけ要素が変位した後、弾性係数などの各ガウスポイント６４の素材特性も変化し得る。各要素５８の、増加分の変位によって起こる素材特性の変化を説明するために、各ガウスポイント６４の弾性係数は変位後アップデートされる。この弾性係数を示すアップデートされた数値の大きさは、グラフ７０の点線部分の上にかぶさる放物線状の部分に示されるように、グラフ的に放物線状の面として図示してある。グラフ７０の放物線状の面はまた、各ガウスポイント６４での剛性への寄与をも図示していて、この寄与は上述の剛性マトリクス［ｋ］を評価するために数学的にモデリングされる。

メッシュが穴のような空の空間を囲む素材上に当てはめられている実施例においては、弾性係数の数値はその穴をモデリングするために故意に小さくしてもよい。例えば、穴を有する金属薄板が素材の例であるならば、穴から遠いガウスポイント６４の弾性係数は２９Ｘ１０^６ＰＳＩと決め、穴に近いガウスポイント６４の弾性係数は１０ＰＳＩと決めることができる。

上述のように弾性係数は、σ＝Ｅεとなるような応力と歪みの間の関係を表す。Ｅは弾性係数を意味し、これは弾性ドメインでの応力と歪みの比率を表す数値である。σは応力を表し、εは歪みを表す。σ＝Ｅεは、歪みの比例限度（ε_ＰＬ）以下の歪みの値に対する有効条件である。要素のドメイン全体の弾性係数は、その例をグラフ７０に図示するが、放物線近似を用いて評価することができ、

であり、またここで、Ｎ_ｉ、Ｅ_ｉ、Ｕ、Ｖ、Ｅ_{（Ｕ、Ｖ）}は、
Ｎ_ｉ＝ノードｉに関連するシェープ関数、
Ｅ_ｉ＝ノードｉにおける弾性係数、
Ｕ＝要素のＵ方向におけるパラメータ値、
Ｖ＝要素のＶ方向におけるパラメータ値、及び
Ｅ_{（Ｕ、Ｖ）}＝要素のパラメータとしてのＵＶポジションにおける弾性係数の補間値である。

要素の各ガウスポイントに対する弾性係数値を、ＦＥＡ法のそれぞれのイテレーションに際し調整することは有益である。何故なら、特定のメッシュ内の要素数を大幅に増加させることなく、デバイスの正確なモデルが生成できるからである。特定のメッシュ内で比較的少ない数の要素を使うと、要求される演算能力と演算時間が少なくてすむ。

図５Ｃは、グローバル変位マトリクス｛Ｄ｝を使いローカル要素座標を確立させるステップ２５８の一つの実施例の更なる詳細を図示する。ステップ３００においては、要素１３０Ａ及び／または１９０Ａに対してシェープ関数（式５乃至６３）を用いた相対位置の初期予測は、これらの要素１３０Ａ及び／または１９０ＡのパラメータとしてのＵＶスペース（パラメータとしてのＵＶスペースにおけるパラメータの例としては−１．０＜＝Ｕ＜＝１．０及び−１．０＜＝Ｖ＜＝１．０がある）に対してのグローバル変位マトリクス｛Ｄ｝と同じで、課された境界条件に基づき、

として評価される。ここで、Ｓ、ｘ、ｙ、ｕ、ｖは、次の通りである。
Ｓ＝要素のｕｖパラメータにおける位置を定義するベクトル。
ｘ＝位置ベクトルＳのＸ座標。
ｙ＝位置ベクトルＳのＹ座標。
ｚ＝位置ベクトルＳのＺ座標。
ｕ＝Ｕ方向での要素のパラメータ値。
ｖ＝Ｖ方向での要素のパラメータ値。

境界条件は、制限された各自由度に対する「ペナルティー法」を用いて課される。ステップ３０４において、直交ローカル座標システムを決定するために、各要素１３０Ａまたは１９０Ａに対する一次空間導関数を評価する。各要素１３０Ａまたは１９０Ａに対する一次空間導関数は、次式として評価される。

ステップ３０４で決定された一次導関数を用い、与えられたＵＶパラメータに対して単位ベクトルe_１、e_２、及びe_３を有する直交座標システムを、次式のようにステップ３０６において計算する。

ここで、ｅ、ｌ、m、nは、次の通りである。
ｅ＝ローカル座標システムの座標ベクトル。
ｌ＝座標ベクトルｅのローカルｘ成分。
m＝座標ベクトルｅのローカルｙ成分。
n＝座標ベクトルｅのローカルｚ成分。
結果として得られるグローバル変位マトリクス｛Ｄ｝は

と表現される。

図５Ｄは、ローカル座標内の要素剛性マトリクスを表す剛性マトリクス［ｋ］及び、ローカル座標の要素力を表す力ベクトル｛ｆ｝を確立／調整するステップ２６４の一つの実施例の更なる詳細を図示する。力ベクトル｛ｆ｝は、計算された剛性マトリクス［ｋ］及びローカルノード変位に基づき決定される。剛性マトリクス［ｋ］は、歪み―変位マトリクス、弾性応力―変位マトリクス及び弾性―塑性応力―変位マトリクス、ヤコビアンマトリクスを用いた、数値積分によって決定される。歪み―変位マトリクスは、各要素に対する変位と歪みの関係を定義する。応力―変位マトリクスは、要素の歪みと応力の関係を定義する。ヤコビアンマトリクスの決定は、各要素の数値積分に必要である。この手順の詳細を、図５Ｄを参照して以下に説明する。

ステップ３３０で、要素１３０Ａ及び１９０Ａに対する応力変換マトリクス及び歪み変換マトリクスが、それぞれ

として定義される。ここで、Ｔ_σ、Ｔ_ε、Ｔ_１１、Ｔ_１２、Ｔ_２１、Ｔ_２２は、
Ｔ_σ＝応力変換マトリクス、
Ｔ_ε＝歪み変換マトリクス、
Ｔ_１１＝左上部象限部分マトリクス、
Ｔ_１２＝右上部象限部分マトリクス、
Ｔ_２１＝左下部象限部分マトリクス、
Ｔ_２２＝右下部象限部分マトリクスであり、
各部分マトリクスが、次式として定義される。

これらの各応力変換マトリクス及び歪み変換マトリクスに対して、

のように、その逆マトリクスは転換マトリクスに等しい。
ステップ３３４において、シェル要素１３０Ａ及び１９０Ａに対してのヤコビアンマトリクスは、以下のように評価される。

ここで、x_ｉ ^＊、ｙ_ｉ ^＊、ｚ_ｉ ^＊は、
ｘ_ｉ ^＊＝要素の垂直方向に沿ったX成分の大きさ、
ｙ_ｉ ^＊＝要素の垂直方向に沿ったＹ成分の大きさ、及び
ｚ_ｉ ^＊＝要素の垂直方向に沿ったZ成分の大きさであり、
またここで、Ｖ_ｉ ^＊が

として定義され、またここで、ｗ、ｔ_ｉは、次の通りである。
ｗ＝要素の厚さ方向に沿ったパラメータの値（−１．０＜＝ｗ＜＝１．０）。
ｔ_ｉ＝要素の厚さ。
そして、然るに逆ヤコビアンが、

として定義される。ここで
ｊ^＊＝ヤコビアンマトリクスの逆要素である。

ステップ３３８において、ε_ｚ’＝０と仮定すると、要素１３０Ａ及び１９０Ａのローカル軸に沿った歪みの成分は以下のように与えられる。

ここで、ε_ｘ’、ε_ｙ’、γ_ｘ’ｙ’、γ_ｘ’ｚ’、γ_ｙ’ｚ’は、次の通りである。
ε_ｘ’＝主ｘ方向の歪み。
ε_ｙ’＝主ｙ方向の歪み。
γ_ｘ’ｙ’＝主ｘｙ面のせん断。
γ_ｘ’ｚ’＝主ｘｚ面のせん断（曲げ）。
γ_ｙ’ｚ’＝主ｙｚ面のせん断（曲げ）。
上記の式は、表面に沿った成分と曲げ成分に関し、

と書き表すことができる。ここで、｛ε_ｍ｝、｛ε_ｓ｝は、

であり、またここで、ε_ｍ、ε_ｓは、
ε_ｍ＝歪みの表面に沿った成分、及び
ε_ｓ＝歪みの曲げ成分である。

ステップ３４０において、各要素に対する歪み変位マトリクスは、以下の数学的プロセスを用いて決定される。グローバル座標において、要素の歪みを計算するために下記の関係を確立することができる。

歪み変位マトリクス［Ｂ］は［Ｂ_１Ｂ_２ΛＢ_ｎ］として与えられるが、ここでＢ_iは、

として定義され、またここで、ａ_ｉ、ｂ_ｉ、ｃ_ｉ、ｄ_ｉ、ｅ_ｉ、ｆ_ｉは、次の通りである。

ローカル座標では、この歪み―変位マトリクスは、

と書ける。

ステップ３４４において、弾性変形時の要素の応力とノード変位は、

によって関連付けることができるが、ここで［Ｃ］は、次の通りである。
［Ｃ］＝構成マトリクス。
σ_ｚ’＝０であると仮定して、ｘ’座標ｙ’座標ｚ’座標における応力―歪み関係に対して次のことが得られる。ｘ’、ｙ’、ｚ’は要素のローカル座標方向を表す。

つまり、

であるが、ここでσ_ｘ’、σ_ｙ’、τ_ｘ’ｙ’、τ_ｙ’ｚ’、τ_ｚ’ｘ’、μは、
σ_ｘ’＝主x方向の応力、
σ_ｙ’＝主y方向の応力、
τ_ｘ’ｙ’＝主ｘｙ面のせん断応力、
τ_ｙ’ｚ’＝主ｚｙ面の曲げ応力、
τ_ｚ’ｘ’＝主ｚｙ面の曲げ応力、及び
μ＝ポアソン率である。

変形可能な物体に軸上の伸張力が加えられると、その物体は伸張するのみならず横方向に収縮もする。同様に、物体に働く軸上の圧縮力は物体がその力の方向に収縮し、その物体の側面は横に広がる。これらの歪みの比例比率はポアソン率（μ）と呼ばれ、ここで用いられる等方性素材では、次の通りとなる。

しかしながら、要素の接線方向に適用されるポアソン率の値と、要素に垂直な方向に適用される値との間の異方性関係を考慮する必要があり得る。この垂直な異方性比率、つまりｒ要因は、

として適用され、ここで、μ_Ｔ、μ_Ｎは、次の通りである。
μ_Ｔ＝表面に沿った方向でのポアソン率。
μ_Ｎ＝厚さ方向でのポアソン率。

等方性のｒ値（ｒ＝１．０）に対しては、実際の体積保存条件を考慮し、ポアソン率の値は、次式の通りとなる。

このようにして、ポアソン率の与えられた等方性値に対して、下記の関係が得られ、下記の関係を用いることで、与えられたｒ値に基づいて対応する異方性値を計算することができる。

この二次方程式を解きμ_Ｎ及びμ_Ｔを求めると、その結果は次の通りとなる。

このようにして、各要素の構成マトリクスを定義する場合は、接線方向のポアソン率の値μTが用いられる。

式８９及び式９４を用い、応力―変位マトリクスは表面に沿った成分及び曲げ成分に分けることができ、

となるが、ここで、Ｃ_m、Ｃ_sは、

であり、またここで、Ｃ_m、Ｃ_sは、
Ｃ_m＝応力―歪み関係の表面に沿った成分の構成マトリクス、及び
Ｃ_s＝応力―歪み関係の曲げ成分の構成マトリクスである。

となるが、ここでＰは歪みの塑性成分を示す。

であり、また、σ_ｘ ^’、σ_ｙ ^’は、次の通りである。

プラントルーロイス関係によると、塑性歪みの増加分は次の式によって定義される。

歪みの増加分（ｄε）は、その弾性成分（ｄε^ｅ）と塑性成分（ｄε^ｐ）を加算したものであるので、次の式の通りである。

となるが、ここでTはベクトルの転置マトリクスを示す。

と書き表すことができる。ここで、［Ｗ］は次の通りである。
［Ｗ］＝歪みの総増加分に関する塑性歪みの増加分を計算するテンソル。

ローカル座標ｘ’、ｙ’、ｚ’での［Ｂ］マトリクスと［ＣＢ］マトリクスを得た後、これらのマトリクスを利用してステップ３５０で要素剛性マトリクスの計算を行う。このマトリクスの変換は、仮想的変位の際に結果として得られる歪みエネルギー密度の増加分がその計算が行われる座標システムに関わらず同じであるべきであるという概念に基づき、導き出すことができる。こうして剛性マトリクスは、

で与えられる。ここで、Ｂ^Ｔ＝歪み―変位マトリクスの転置マトリクスである。

一つの実施例においては、四角形シェル要素１３０Ａに対して、下記のガウスポイントと重さを用いた３ｘ３ｘ２の次数のガウス―ルジャンドル求積法を使って積分を行うことができる。

ここで、Ｗ＝ガウスポイントにおけるウェイト要因である。

一つの実施例においては、三角形シェル要素１９０Ａに対して、下記のガウスポイントと重さを用いて、積分を行ことができる。

その後ステップ３５２において、力ベクトル｛ｆ｝が式｛ｆ｝＝−［Ｋ］｛ｄ｝を用いて、剛性マトリクス［Ｋ］に基づき決定される。

一つの実施例においては、最初に領域５４の塑性変形解析を行った後に、方法２５０のステップ２８８でオプションとしてスプリングバック解析及び／またはオーバークラウン解析を行うことができる。スプリングバック解析とは、それまでに適用された境界条件が解除され素材がリラックスした、つまり自然な状態に戻った後の素材の形状の解析を意味する。オーバークラウン解析とは、スプリングバックした形状が素材の望ましい形状となるために必要な、素材の過剰成形の度合いを予測する事を意味する。スプリングバック／オーバークラウン解析は、素材の形成過程のために生じる内部応力の蓄積を説明する。スプリングバック／オーバークラウン解析は、未形成部分から組み立てられる剛性マトリクスに基づく静的解析を用いて行うことができる。モデルに加えられるノードでの力は、方法２５０において記述した最初の形成過程の際に得られた修正されたグローバル剛性マトリクスの逆マトリクスから決定され、これは

及び

である。ここで、｛Ｆ｝、［Ｋ_Ｉ］^−１、［Ｋ_Ｆ］、｛Ｄ_Ｐ｝、｛Ｄ_Ｓ｝は次の通りである。
｛Ｆ｝＝修正されたグローバル剛性マトリクスの逆マトリクスと、初期形成解析の結果として得られる変位とから計算された、結果として得られる内部負荷ベクトル。
［Ｋ_Ｉ］^−１＝修正されたグローバル剛性マトリクスの逆マトリクス。
［Ｋ_Ｆ］＝初期形成解析から結果として得られた変位したノード座標から計算された剛性マトリクス。
｛Ｄ_Ｐ｝＝初期形成解析から計算された、結果として得られた変位ベクトル。
｛Ｄ_Ｓ｝＝ポスト線形静的解析から計算されたスプリングバック変位ベクトル。
このアプローチは、形状形成された素材の弾性―スプリングバック形状を誘導するために必要な応力の正反対のものに基づくノードの負荷［Ｋ_Ｉ］^−１を適用することによって、オーバークラウンの形状は合理的に近似されている、と仮定する。

本発明は詳細に説明されているが、頭記の請求項の範囲に定義される本発明の精神及び範囲から逸脱せずに、様々な変更、置き換え、または代替が行うことが可能であることが理解されるべきである。

添付した図に併せた以下の説明についてここで言及する。添付の図においては、同一の参照番号は同一のパーツを示すものとする。

有限要素解析の一つの方法で使用できる初期のメッシュの実施例を図示する概略図である。 目的の形状を有し、有限要素解析の方法を実行した結果得られる最終メッシュの実施例を図示する概略図である。 モデリングされる例としての素材に関連する応力と歪みの関係のグラフの例である。 本発明の一つの実施例による有限要素解析の一つの方法を実行するために作動可能なコンピュータシステムの実施例を図示する概略図である。 図３Ａに示すコンピュータの実施例を図示するブロック構成図である。 図２Ａに示す初期のメッシュの要素の実施例の概略図である。 図２Ａに示す初期のメッシュの要素の実施例の概略図である。 本発明の一つの実施例による有限要素解析の方法を図示するフローチャートである。 図５Ａの方法を実施するために使用できる弾性係数の放物線近似を図示するグラフである。 図５Ａの方法に示されるステップの更なる詳細を図示するフローチャートである。 図５Ａの方法に示されるステップの更なる詳細を図示するフローチャートである。

Claims (28)

1. 金属薄板の表面の変形を決定するコンピュータ化された方法であって、
   コンピュータが、
   該金属薄板の弾性変形範囲及び塑性変形範囲を受信するステップと、
   複数の要素から成るメッシュで該金属薄板の一領域をモデリングするステップと、
   コンピュータが、
   該メッシュに対する有限要素解析を行うことによって、要素のそれぞれ上の少なくとも一つのポイントに対する変位を決定するステップと、
   なお、該変位はモデリングされた一連の境界条件の結果として得られ、
   コンピュータが、
   以前に計算された要素の変位から計算される歪みの値に基づいて該弾性変形領域及び該塑性変形領域に従って各要素に対する弾性係数を反復的にアップデートするステップと、を有し、
   更に、
   前記金属薄板の弾性変形範囲と該塑性変形範囲を受信するステップが、降伏応力に関連する値を受信するステップであり、
   前記各要素に対する該弾性係数を反復的にアップデートするステップが、
   以前に決められた要素の変位から計算される歪みの値が前記降伏応力に関連する値よりも小さい場合には、前記反復的にアップデータされた弾性係数が、該弾性変形範囲に従って弾性係数とするステップであり、
   コンピュータが、
   該アップデートされた弾性係数を使い該剛性マトリクスを計算するステップと、
   該剛性マトリクスは、歪み―変位マトリクス、弾性応力―変位マトリクス及び弾性―塑性応力―変位マトリクス、ヤコビアンマトリクスを用いた、数値積分によって決定される、を有することを特徴とする方法。
2. 該複数の要素の該それぞれが少なくとも一つのノードから成ることを特徴とする請求項１に記載の方法。
3. 前記方法が、
   コンピュータが、該複数の要素の該それぞれに対して、剛性マトリクス及び少なくとも一つの力ベクトルを該弾性係数に基づき決定するステップを更に有する、
   なお、該剛性マトリクス及び該少なくとも一つの力ベクトルが、該複数の要素の該それぞれに特有なローカル座標系のフォーマットで表される、請求項１に記載の方法。
4. 前記方法が、
   コンピュータによって、
   複数の要素の該それぞれに対して、剛性マトリクス及び少なくとも一つの力ベクトルを弾性係数に基づき決定されるステップと、
   なお、該剛性マトリクス及び該少なくとも一つの力ベクトルが、該複数の要素の該それぞれに特有なローカル座標系のフォーマットで表され、
   コンピュータによって、
   該有限要素解析を行うことが、該剛性マトリクス及び該少なくとも一つの力ベクトルを、該複数の要素の全てに対して使用されるグローバル座標系のフォーマットへ変換されるステップと、
   ことを更に含むことと、を特徴とする請求項１に記載の方法。
5. その時点での変位と以前の変位との差が事前に決められたモデリングの許容範囲よりも小さくなった後に、
   コンピュータが、
   該複数の要素のそれぞれに対応する複数の剛性マトリクスの合計から、グローバル剛性マトリクスを決定するステップと、
   コンピュータが、
   該複数の要素の該それぞれに対して、
   ｛Ｄｓ｝＝［Ｋｆ］・［Ｋｌ］^−１・｛Ｄｐ｝
   を計算することによって、スプリングバックのポジション｛Ｄｓ｝を決定するステップとを有し、
   ただし、［Ｋ_Ｉ］^−１＝修正されたグローバル剛性マトリクスの逆マトリクス、
   ［Ｋ_Ｆ］＝初期形成解析から結果として得られた変位したノード座標から計算された剛性マトリクス、
   ｛Ｄ_Ｐ｝＝初期形成解析から計算された、結果として得られた変位ベクトル、である、
   ことを特徴とする請求項１に記載の方法。
6. 該複数の要素の少なくとも一つが四角形であり、少なくとも１２個のノードを該四角形のノードのラインに沿って指定することから更に成ることと、該少なくとも１２個のノードの内の４個がそれぞれ該四角形の４個の角に配置されていることと、を特徴とする請求項１に記載の方法。
7. コンピュータが
   該金属薄板の降伏応力に関連するポイントにおいて該塑性変形範囲に変換される該弾性変形範囲を受信するステップと、
   各要素に対する該弾性係数を反復的にアップデートする場合であって、該以前に計算された要素の変位から計算される歪みの該値が、降伏応力のポイントに達している場合、
   コンピュータによって、
   前記アップデートされる弾性係数が、その反応として、該弾性係数が該塑性変形範囲内にある弾性係数としてアップデートされるステップを更に有する、請求項１に記載の方法。
8. 該弾性変形範囲と該塑性変形範囲を受信することが降伏応力に関連する値を受信することを含むことと、
   各要素に対する該弾性係数を反復的にアップデートすることが、以前に計算された要素の歪曲から計算された歪みの該値が降伏応力に関連する値に達していなと決定し、その反応として、該弾性係数が該弾性変形範囲内にあると決定することから成ることと、
   を特徴とする請求項１に記載の方法。
9. 該弾性係数をアップデートすることが、該要素に課された実際の応力を利用することから成り、該実際の応力は、フォンミーゼスの条件を使って計算されることを特徴とする請求項１に記載の方法。
10. 有限要素解析の方法を使用して、金属薄板の表面の変形を決定するコンピュータ化された方法であって、
    コンピュータが、
    該金属薄板の弾性変形範囲と塑性変形範囲を受信するステップと、
    複数の要素から成るメッシュで該金属薄板をモデリングするステップと、
    なお、該複数の要素のそれぞれが複数のノードから成り、
    コンピュータが、
    或る時点での変位と、その前のステップでの変位との差が事前に決められたモデリングの許容範囲よりも小さくなるまで繰り返される、有限要素解析の計算ステップの各ステップにおいて、該複数の要素の該それぞれに対して剛性マトリクスを計算するステップと、
    コンピュータが、
    その前のステップで決められた要素の変位から計算される歪みの値が、降伏応力に関連する値よりも小さい場合には、
    前記反復的にアップデータされた、該弾性変形範囲に従って弾性係数を使い該剛性マトリクスを計算するステップと、
    および、
    該その前のステップで決められた要素の変位から計算される歪みの値が、降伏応力に関連する値と等しいか、またはそれよりも大きい場合には、
    前記反復的にアップデータされた、該塑性変形範囲に従って弾性係数を使い、該剛性マトリクスを計算するステップと、
    コンピュータが
    該複数のノードの該それぞれの変位が事前に決められたモデリングの許容範囲内に収まるまで、該複数の要素の該それぞれに対する該剛性マトリクスを計算するステップと、
    該剛性マトリクスは、歪み―変位マトリクス、弾性応力―変位マトリクス及び弾性―塑性応力―変位マトリクス、ヤコビアンマトリクスを用いた、数値積分によって決定される
    上記各ステップを有することを特徴とする方法。
11. 該複数の要素のそれぞれが少なくとも一つのガウスポイントから成ることと、剛性マトリクスを計算することが、アップデートされた弾性係数を複数のガウスポイントのそれぞれに割り当てることから成り、該アップデートされた弾性係数が該複数の要素の該それそれの歪みの値の変化に従って決定され、該値の該変化が、前記繰り返される有限要素解析の、或る時点での変位と、その前のステップでの変位との差である増加分を基に計算されることと、を特徴とする請求項１０に記載の方法。
12. 該弾性変形範囲と該塑性変形範囲のそれぞれを、歪みと応力の直線関係として近似することから更に成ることを特徴とする請求項１０に記載の方法。
13. 該剛性マトリクスの計算が、該剛性マトリクスのフォーマットを該複数の要素の全てに継続的に適用するために作動可能であるグローバル座標系を用いて、新しいフォーマットに変換することから更に成ることを特徴とする請求項１０に記載の方法。
14. 該複数の要素の該それぞれに関連する少なくとも一つの力ベクトルを、該複数の要素のそれぞれに対応する該剛性マトリクスを用いて計算し、該少なくとも一つの力ベクトルと、該複数の要素の該それぞれの該剛性マトリクスとが、それぞれに対応するローカル座標系にあることと、
    該剛性マトリクス及び該要素のそれぞれの該少なくとも一つの力ベクトルを、該複数の要素の全てに対して使用されるグローバル座標系のフォーマットに変換することと、
    から更に成ることを特徴とする請求項１０に記載の方法。
15. その時点での変位と以前の変位との差が事前に決められたモデリングの許容範囲よりも小さくなった後に、
    コンピュータが、
    該複数の要素のそれぞれに対応する複数の剛性マトリクスの合計からグローバル剛性マトリクスを決定するステップと、
    該複数の要素の該それぞれに対して、
    ｛Ｄｓ｝＝［Ｋｆ］・［Ｋｌ］^−１・｛Ｄｐ｝
    を計算することによって、スプリングバックのポジション｛Ｄｓ｝を決定するステップと、を有する、
    ただし、［Ｋ_Ｉ］^−１＝修正されたグローバル剛性マトリクスの逆マトリクス、
    ［Ｋ_Ｆ］＝初期形成解析から結果として得られた変位したノード座標から計算された剛性マトリクス、
    ｛Ｄ_Ｐ｝＝初期形成解析から計算された、結果として得られた変位ベクトル、である、
    ことを特徴とする請求項１０に記載の方法。
16. 該複数の要素の少なくとも一つが四角形であり、該四角形のノードのラインに沿って少なくとも１２個のノードを指定することから更に成ることと、該少なくとも１２個のノードの内の４個がそれぞれ該四角形の４個の角に配置されていることと、を特徴とする請求項１０に記載の方法。
17. 剛性マトリクスを計算するときに、該要素に課された実際の応力に基づいて弾性係数をアップデートする、ここで、該実際の応力は、フォンミーゼスの条件を用いて計算されることを特徴とする請求項１０に記載の方法。
18. 弾性変形範囲が、該塑性変形範囲の応力対歪みの関係とは異なる応力対歪みの関係を有することを特徴とする請求項１０に記載の方法。
19. 金属薄板の変位のモデリングさせるために、コンピュータに下記ステップを実行させるためのプログラムであって、
    コンピュータが、
    該金属薄板の弾性変形範囲及び塑性変形範囲を受信するステップと、
    複数の要素から成るメッシュで該金属薄板の一領域をモデリングするステップと、
    を実行し、
    該コンピュータが、
    該メッシュの有限要素解析を行うことによって、該要素のそれぞれ上の少なくとも一つのポイントに対する変位を決定するステップと、
    コンピュータが、
    以前に計算された要素の変位から計算される歪みの値に基づいて該弾性変形領域及び該塑性変形領域に従って各要素に対する弾性係数を反復的にアップデートするステップと、
    更に、
    前記金属薄板の弾性変形範囲と該塑性変形範囲を受信するステップが、降伏応力に関連する値を受信するステップと
    を実行し、
    前記各要素に対する該弾性係数を反復的にアップデートするステップが、
    以前に決められた要素の変位から計算される歪みの値が前記降伏応力に関連する値よりも小さい場合には、前記反復的にアップデータされた弾性係数が、該弾性変形範囲に従って弾性係数とするステップであり、
    コンピュータが、
    該アップデートされた弾性係数を使い該剛性マトリクスを計算するステップと、
    該剛性マトリクスは、歪み―変位マトリクス、弾性応力―変位マトリクス及び弾性―塑性応力―変位マトリクス、ヤコビアンマトリクスを用いた、数値積分によって決定される、
    を実行する、ことを特徴とするプログラム。
20. 該複数の要素の該それぞれが少なくとも一つのノードから成ることを特徴とする請求項１９に記載のプログラム。
21. 該プログラムが、該複数の要素の該それぞれに対して、剛性マトリクス及び少なくとも一つの力ベクトルを該弾性係数に基づき決定することによって、該有限要素解析を行うために作動可能であることと、該剛性マトリクス及び該少なくとも一つの力ベクトルが、該複数の要素の該それぞれに特有なローカル座標系のフォーマットであることと、を特徴とする請求項１９に記載のプログラム。
22. コンピュータによって、
    該複数の要素の該それぞれに対して、剛性マトリクス及び少なくとも一つの力ベクトルを弾性係数に基づき決定されるステップと、
    なお、該剛性マトリクス及び該少なくとも一つの力ベクトルが、該複数の要素の該それぞれに特有なローカル座標系のフォーマットで表され、
    コンピュータによって、
    該剛性マトリクス及び該少なくとも一つの力ベクトルを、該複数の要素の全てに対して使用されるグローバル座標系のフォーマットへ変換されるステップと、を更に有する、請求項１９に記載のプログラム。
23. その時点での変位と以前の変位との差が事前に決められたモデリングの許容範囲よりも小さくなった後に、
    コンピュータが、
    該複数の要素のそれぞれに対応する複数の剛性マトリクスの合計からグローバル剛性マトリクスを決定するステップと、
    該複数の要素の該それぞれに対して、
    ｛Ｄｓ｝＝［Ｋｆ］・［Ｋｌ］^−１・｛Ｄｐ｝
    を計算することによって、スプリングバックのポジション｛Ｄｓ｝を決定するステップと、を有する、
    ただし、［Ｋ_Ｉ］^−１＝修正されたグローバル剛性マトリクスの逆マトリクス、
    ［Ｋ_Ｆ］＝初期形成解析から結果として得られた変位したノード座標から計算された剛性マトリクス、
    ｛Ｄ_Ｐ｝＝初期形成解析から計算された、結果として得られた変位ベクトル、である、
    請求項１９に記載のプログラム。
24. 該複数の要素の少なくとも一つが四角形であり、少なくとも１２個のノードを該四角形のノードのラインに沿って指定することから更に成ることと、該少なくとも１２個のノードの内の４個がそれぞれ該四角形の４個の角に配置されていることと、を特徴とする請求項１９に記載のプログラム。
25. コンピュータによって、
    該金属薄板の降伏応力に関連するポイントにおいて該弾性変形範囲が該塑性変形範囲に変換されることと、
    コンピュータによって、
    該以前に計算された要素の変位から計算される歪みの該値が降伏応力のポイントに達していると決定されるステップと、
    コンピュータによって、
    それぞれの要素に対して該弾性係数をアップデートされるステップと、
    なお、該弾性係数が該塑性変形範囲内にあると決定されるステップと、更に有する、請求項１９に記載のプログラム。
26. コンピュータによって、
    該弾性変形範囲と該塑性変形範囲との間の変換ポイントである該金属薄板の降伏応力に関連する値を受信されるステップと、
    コンピュータによって、
    以前に計算された要素の歪曲から計算された歪みの該値が降伏応力に関連する値に達していない場合、各要素に対する、該弾性変形範囲内にある、該弾性係数を反復的にアップデートするステップと、を更に有する
    を特徴とする請求項１９に記載のプログラム。
27. コンピュータによって、
    該要素に課された実際の応力を利用することによって該弾性係数をアップデートするステップを有する、
    なお、該実際の応力は、フォンミーゼスの条件を使って計算される、
    ことを特徴とする請求項１９に記載のプログラム。
28. コンピュータ支援解析環境下での金属薄板の変位をモデリングするためのシステムであって、
    該システムが、
    ディスプレイユニット、入力デバイス、及びプロセッサを有するコンピュータシステムと、
    該コンピュータシステムに結合された、コンピュータで読み取り可能な媒体と、から成り、該コンピュータで読み取り可能な媒体は、該プロセッサ上で実行するために作動可能なソフトウェアプログラムから成り、
    前記ソフトウエアプログラムは、前記請求項１９乃至２７いずれか１項記載のプログラムであることを特徴とするシステム。

Images