JPWO2010038539A1 - 成形シミュレーション方法、成形シミュレーション装置、プログラム、記録媒体、及びシミュレーション結果に基づいた成形方法 - Google Patents

Abstract

この弾塑性材料の成形シミュレーション方法は、有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する工程と；計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する工程と；を含む。

Description

本発明は、成形シミュレーション方法、成形シミュレーション装置、プログラム、記録媒体、及びシミュレーション結果に基づいた成形方法に関する。本発明は特に、鋼板やアルミ板等の金属板や、プラスチック材料や複合材等の弾塑性材料を成形する際の成形シミュレーション及びその結果に基づいた成形方法に関する。 本願は、２００８年９月３０日に、日本に出願された特願２００８−２５２７３０号、及び２００８年９月３０日に、日本に出願された特願２００８−２５３２７８号に基づき優先権を主張し、それらの内容をここに援用する。

従来より、自動車部品や家庭電気製品には、上下一対の凹凸金型を用いて、薄鋼板やアルミ薄板等の金属薄板を所定の形状にプレス加工した部品が多用される。近年、このような製品の部品の軽量化を図るため、より高強度の材料を用いることで板厚を減少させることが盛んに行われている。一方、高強度の材料をプレス成形する際には、材料の弾性変形が発生し、成形下死点における部品形状と異なる形状に製品が形成されてしまうことが多い。このスプリングバック（弾性回復）の量が大きくなると、部品精度の低下を招く。これを抑制するために、部品の形状を変更することで部品剛性を高めることにより、スプリングバック量を低減することが通常行われている。また、スプリングバック量を予め見込んで、金型形状を修正することも、通常行われている。

一方、プレス部品は一般に複雑な形状を有しているため、成形時の歪みは一様ではない。従って、スプリングバックの発生原因部位を特定することは簡単ではない。これに対して、対象部品の成形下死点における応力分布、すなわち対象部品の目標形状に基づき、部品を複数の領域に分割し、それぞれの領域における応力を順次変化させて、スプリングバック解析を行い、スプリングバックに支配的な領域を特定する方法が下記特許文献１〜３に開示されている。これら特許文献１〜３に記載の方法は、いずれも、非特許文献１に記載の有限要素法を利用した方法である。

さらに、スプリングバックの原因となる内部残留応力を制御する方法として、例えば特許文献４〜７には、特徴的な部品形状に対して、予め特定の部位にエンボスや突起形状を成形し、次工程で潰す方法が開示されている。また特許文献８には、素材（ブランク）全面にエンボスや突起形状を成形し、次工程で潰す方法が開示されている。

特開２００７−２２９７２４号公報 特開２００８−４９３８９号公報 特開２００８−５５４７６号公報 特開２００６−２７２４１３号公報 特開２００７−２２２９０６号公報 特開２００８−１２５７０号公報 特開２００８−１８４４２号公報 特開２００６−３５２４５号公報

有限要素法ハンドブック：鷲津久一郎 他編、培風館（１９８１）

しかしながら、特許文献１〜３の手法では、いずれも、予め部品を複数の領域に分割する。そして、分割された領域の数だけ、大規模な連立方程式求解であるスプリングバック計算を繰り返し行う必要がある。従って、スプリングバック発生の原因となる部位を特定する作業が煩雑化する。また、分割の仕方（大きさ、分割数）により、結果が異なる。従って、スプリングバック発生の原因となる部位を十分に特定することが困難であるという問題がある。

本発明は、上記の課題に鑑みてなされたものであり、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速かつ確実に、金属板の成形時におけるスプリングバックの発生部位を効率良く特定し、正確な成形に供することを可能とする成形シミュレーション方法、成形シミュレーション装置、プログラム、並びに記録媒体を提供することを第１の目的とする。

また、特許文献４〜７の方法は、いずれも、Ｕの字型やハット型断面で、高さ方向あるいは長手方向に湾曲する比較的単純形状部材の、ウェブ面やフランジ面にエンボスや突起形状を成形し、次工程で潰す。これに対し、実際の例えば自動車部品などでは、一般に形状は複雑で、ウェブ面やフランジ面に他部品との接合などの目的で凹凸形状や開口部がある。従って成形時の応力や歪みはウェブ面内やフランジ面内で一様でなく、複雑に分布している。このため、特許文献４〜７の方法では、エンボスを付与する部位によっては効果が得られない、あるいはスプリングバックを増大させてしまう、といった問題がある。

また特許文献８の方法では、ブランク全面にエンボスを付与するため、潰すための成形力が増大するという問題がある。更に、エンボスを潰すことにより圧縮応力が発生し、これが駆動力となってスプリングバックを増大させる、といった問題もある。
また、特許文献４〜７の方法と特許文献１に記載の方法とを組み合わせた場合であっても、複雑形状のプレス部品のプレス加工時におけるスプリングバックを十分に抑制することはできなかった。

本発明は、上記の課題に鑑みてなされたものであり、複雑形状のプレス部品に対して、プレス成形時におけるスプリングバックの発生の原因となる部位を効率良く特定し、その部位に起因するスプリングバックを抑えることで寸法精度に優れた成形方法を提供することを第２の目的とする。

上記の目的を達成するために、本発明は以下の手段を採用した。
（１）本発明の弾塑性材料の成形シミュレーション方法は、有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する工程と；計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する工程と；を含む。
（２）前記（１）の成形シミュレーション方法は、計算された前記全等価節点力ベクトルに基づいて、前記弾塑性材料の全領域又は特定の領域のうちで、前記等価節点力ベクトルの大きい部位をスプリングバックの発生原因部位として特定する工程を更に含んでも良い。
（３）前記（１）に記載の成形シミュレーション方法は、全体剛性行列の逆行列を計算する工程と；前記全等価節点力ベクトルを外力ベクトルとして、前記全体剛性行列の逆行列と前記全等価節点力ベクトルとを乗算し、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの寄与度を計算する工程と；を更に含んでも良い。
（４）前記（３）に記載のシミュレーション方法は、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの前記寄与度を表示する工程を更に含んでも良い。
（５）前記（１）に記載の成形シミュレーション方法は、全体剛性行列の逆行列を計算する工程と；前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する工程と；前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する工程と；前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する工程と；を更に含んでも良い。
（６）前記（５）に記載の成形シミュレーション方法は、前記第１の変位量と前記第２の変位量との間の変化量を表示する工程を更に含んでも良い。
（７）本発明の弾塑性材料の成形方法は、前記（１）に記載のシミュレーション方法の結果に基づいて、弾塑性材料を成形する。
（８）前記（７）に記載の弾塑性材料の成形方法は、前記シミュレーション方法に基づいてスプリングバック発生原因部位を特定する工程と；前記スプリングバック発生原因部位に予め、エンボスを成形する工程と；前記エンボスに対し、圧縮応力を与えるように塑性変形させる工程と；を含んでも良い。
（９）前記（８）に記載の弾塑性材料の成形方法は、前記エンボス成形工程で成形する前記エンボスが、２個以上の略同一寸法の円形エンボスであっても良い。
（１０）前記（９）に記載の弾塑性材料の成形方法は、前記エンボス成形工程が、金型に設けられた凹部に勘合する高さ寸法の異なる複数の着脱式の工具を用いてエンボス形状の高さ寸法を調整する工程を有しても良い。
（１１）本発明のプログラムは、コンピュータにより、弾塑性材料の成形シミュレーションを行うに際して、有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する工程と；計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する工程と；をコンピュータに実行させる。
（１２）前記（１１）に記載のプログラムは、計算された前記全等価節点力ベクトルに基づいて、前記弾塑性材料の全領域又は特定の領域のうちで、前記等価節点力ベクトルの大きい部位をスプリングバックの発生原因部位として特定する工程を更にコンピュータに実行させても良い。
（１３）前記（１２）に記載のプログラムは、全体剛性行列の逆行列を計算する工程と；前記全等価節点力ベクトルを外力ベクトルとして、前記全体剛性行列の逆行列と前記全等価節点力ベクトルとを乗算し、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの寄与度を計算する工程と；を更にコンピュータに実行させても良い。
（１４）前記（１１）に記載のプログラムは、全体剛性行列の逆行列を計算する工程と；前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する工程と；前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する工程と；前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する工程と；を更にコンピュータに実行させても良い。
（１５）本発明のコンピュータで読み取り可能な記録媒体は、前記（１１）に記載のプログラムを記録している。
（１６）本発明の、弾塑性材料の成形シミュレーションを行う成形シミュレーション装置は、有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する第１の計算部と；計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する第２の計算部と；全体剛性行列の逆行列を計算する第３計算部と；前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する第４の計算部と；前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する第５の計算部と；前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する第６の計算部と；を備える。

前記（１）乃至（２）の発明によれば、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速で確実に、弾塑性材料の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。
前記（３）乃至（４）の発明によれば、前記（１）乃至（２）の発明よりも更に確実に弾塑性材料の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。
前記（５）乃至（６）の発明によれば、前記（３）乃至（４）の発明よりも更に確実に弾塑性材料の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。
前記（７）乃至（１０）の発明によれば、（１）乃至（６）に記載のシミュレーション方法の結果に基づいて、複雑形状の加工部品であってもプレス成形時におけるスプリングバックの発生部位を効率良くピンポイントに特定し、その部位に起因するスプリングバックを抑えることで、寸法精度に優れた弾塑性材料の成形方法を提供することが可能となる。

本実施例の成形シミュレーションの対象となる金属薄板の成形目標形状（成形下死点における形状）の一例を示す概略斜視図である。 ４節点シェル要素を示す模式図である。 第１実施形態による成形シミュレーション装置の概略構成を示すブロック図である。 第１実施形態による成形シミュレーション方法をステップ順に示すフロー図である。 第１実施形態に用いる金属部品（高張力鋼板）を示す模式図である。 図５Ａの、Ａで示す箇所の部分拡大図である。 成形シミュレーションによって得られた金属部品の離型後のスプリングバック変位量の分布を示す模式図である。 第１実施形態において、表示部による画像表示の一例を示す模式図である。 第１実施形態において、表示部による画像表示の一例を示す模式図である。 第１実施形態の比較例による金属部品の変形結果を示す模式図である。 第２実施形態による成形シミュレーション装置の概略構成を示すブロック図である。 第２実施形態による成形シミュレーション方法をステップ順に示すフロー図である。 第２実施形態に用いる金属部品（高張力鋼板）を示す模式図である。 図１１ＡのＡ部拡大図である。 成形シミュレーションによって得られた金属部品の離型後のスプリングバック変位量の分布を示す模式図である。 第２実施形態において、表示部による画像表示の一例を示す模式図である。 第２実施形態において、表示部による画像表示の一例を示す模式図である。 第２実施形態の比較例による金属部品の変形結果を示す模式図である。 第３実施形態による成形シミュレーション装置の概略構成を示すブロック図である。 第３実施形態による成形シミュレーション方法をステップ順に示すフロー図である。 第３実施形態に用いる金属部品（高張力鋼板）を示す斜視図である。 図１７Ａの部分拡大図である。 成形シミュレーションによって得られた金属部品の離型後のスプリングバック変位量の分布を示す模式図である。 第３実施形態において、表示部による画像表示の一例を示す模式図である。 第３実施形態の比較例による金属部品の変形結果を示す模式図である。 第４実施形態における、パーソナルユーザ端末装置の内部構成を示す模式図である。 プレス成形の対象となる金属薄板プレス部品の一例を示す概略斜視図である。 成形シミュレーションによって得られた金属部品の離型後のスプリングバック変位量の分布を示す模式図である。 成形シミュレーションによって得られた金属部品の成形下死点における応力分布を示す模式図である。 第５実施形態による成形シミュレーションによって得られた金属部品の捩れに対する影響度の分布を示す模式図である。 第５実施形態によるプレス第一工程成形部品の平面図である。 第５実施形態によるプレス第一工程成形部品および金型の図２６に対するＡ−Ａ断面図である。 第５実施形態によるプレス第二工程成形部品の平面図である。 第５実施形態によるプレス第二工程成形部品および金型の図２８に対するＡ−Ａ断面図である。 第５実施形態によるエンボス成形用円筒工具を示す模式図である。

以下、金属薄板のプレス成形シミュレーションを例にとって本発明の実施形態を説明するが、本発明の適応対象は金属薄板をプレス成形するシミュレーションのみに限られるものではない。本発明は、プラスチック材や複合材料等の弾塑性材料を成形する場合や、ロール成形等により成形する場合のシミュレーション等にも適用できる。また、薄板材料の成形のみに限らず、線材や、ある程度厚みのある材料の成形シミュレーション等にも適用できる。

以下に、本発明の第１実施形態について、図１〜８を参照しながら詳細に説明する。

図１は、本発明の成形シミュレーションの対象となる金属薄板の成形下死点における形状（成形目標形状）の一例を示す概略斜視図である。

図１では、金属薄板１０の一部に有限要素のメッシュＭを示している。有限要素としては、ここでは図２に示すような４節点シェル要素を用いる。この４節点シェル要素は、以下の式（１）のように各節点が全体座標系で参照される６つの自由度を有する。また、各積分点は平面応力の以下の式（２）のように３つの成分を有する。また図示されていないが、金属薄板の厚み方向には数層（ここでは５層）の積分点が存在する。

先ず、各有限要素に対する要素等価節点力ベクトルを、以下の式（３）を用いて計算する。ここでは１つの有限要素ごと要素等価節点力ベクトルを計算するが、複数の有限要素からなる有限要素グループごとに前記要素等価節点力ベクトルを計算しても良い。

ここで、［Ｌ］は座標変換マトリックスである。また、［Ｂ］は変位−歪み関係マトリックスである。座標変換マトリックス［Ｌ］は、要素座標系での要素等価節点力ベクトルを全体座標系に変換する変換行列である。この座標変換マトリックス［Ｌ］は、全体座標系の各Ｘ，Ｙ，Ｚ軸の要素座標系（Ｘ‘，Ｙ’，Ｚ‘）に対する余弦を用いて計算される。また、変位−歪み関係マトリックス［Ｂ］は、用いる有限要素のタイプによって異なり、例えば非特許文献１に各種有限要素に対する定式化が示されている。

続いて、以下の式（４）に示すように、要素毎に計算された要素等価節点力ベクトル｛ｆ｝_ｅを、金属薄板の全領域に対して積分し、全等価節点力ベクトル｛ｆ｝を計算する。ここでは金属薄板の全領域に対して積分するが、金属薄板の特定領域（例えば、特に形状が複雑な領域）に対して積分し、その特定領域の全等価節点力ベクトルを計算しても良い。

算出された全等価節点力ベクトルは、成形下死点における応力を等価な節点力（内力）に変換したものとなっており、全体座標系における、断面力（Ｎｘ，Ｎｙ，Ｎｚ）及び曲げモーメント（Ｍｘ，Ｍｙ，Ｍｚ）からなる。
本発明者らは、多数の部品事例で評価を行った。その結果、本発明者らは、等価節点力のうち、曲げモーメント３成分がスプリングバックの駆動力となることを発見した。更に、本発明者らは、応力分布とは異なり、曲げモーメント量の大きい領域は、金属薄膜の限定的な一部分に局所化することを発見した。

上記の発見に基づき、本発明の第１実施形態では、算出された全等価節点力ベクトルの曲げモーメントについて、等高線図等を用いて曲げモーメント量の大きな領域を表示する。この構成により、例えば大規模な連立方程式の求解のような煩雑な作業を一度も行うことなく、極めて短時間にスプリングバック発生の原因となる部位を言わばピンポイントに正確に特定することができる。

（第１実施形態による成形シミュレーション装置）
図３は、本第１実施形態による成形シミュレーション装置の概略構成を示すブロック図である。

この成形シミュレーション装置は、有限要素法を用いて金属薄板の成形下死点における成形シミュレーションを行うものであって、各要素等価節点力ベクトルを算出する第１の計算部１０１と、全等価節点力ベクトルを算出する第２の計算部１０２と、金属薄板におけるスプリングバックの発生原因部位を特定する特定部１０３と、特定部１０３の特定結果を表示する表示部１０４とを有して構成されている。ここで、第１の計算部１０１、第２の計算部１０２、及び特定部１０３は、例えばコンピュータの中央処理装置（ＣＰＵ）の各機能として実現される。

第１の計算部１０１は、与えられた応力テンソル（式（２）参照）に基づいて、例えば式（３）の演算を行い、各有限要素（例えば４節点シェル要素）ごとに要素等価節点力ベクトルを計算する。尚、複数の有限要素を１グループとして、グループごとに要素等価節点ベクトルを計算しても良い。
第２の計算部１０２は、第１の計算部１０１により計算された有限要素ごとの要素等価節点力ベクトルを金属薄板の全領域に対して積分し（式（４）参照）、全等価節点力ベクトルを計算する。尚、ここでは全領域に対して要素等価節点力ベクトルを積分したが、所定の領域、例えば形状が複雑な領域に対して要素等価節点力ベクトルを積分しても良い。

特定部１０３は、第２の計算部１０２により計算された全等価節点力ベクトルの成分のうち、曲げモーメント（Ｍｘ，Ｍｙ，Ｍｚ）について、金属薄板の成形下死点における領域のうちで、曲げモーメント量の大きい部位の位置を、金属薄板におけるスプリングバックの発生原因部位として特定する。
曲げモーメント量の大きい部位の位置の具体的な特定方法としては、例えば、予め規定された所定の閾値と各領域の曲げモーメント量とを比較する。次に、その閾値よりも大きい曲げモーメント量を示す領域を、金属薄板におけるスプリングバックの発生原因部位として特定する。ここで、相異なる複数の閾値を規定しておき、複数の閾値のうちから適宜な閾値を用いるようにしても良い。

表示部１０４は、曲げモーメント量を、例えば金属薄板の形状に対応させて画像として表示する。具体的には、例えば金属薄板の全体像をモーメント量に応じて色付けし、モーメント量の大きいスプリングバックの発生原因部位を視認できるように表示する。
ここで、表示部１０４を、上記の機能に代えて、或いは上記の機能と共に、特定部１０３の特定結果を金属薄板の形状に対応させて画像表示する機能を有するように構成しても好適である。

（第１実施形態による成形シミュレーション方法）
以下、上記の成形シミュレーション装置を用いた成形シミュレーション方法について説明する。

図４は、第１実施形態による成形シミュレーション方法をステップ順に示すフロー図である。
金属薄板としては、ここでは図５Ａ及び図５Ｂに示す金属部品（高張力鋼板）を成形シミュレーションの対象とする。ここで、図５Ａが金属部品の全体図、図５Ｂが部分拡大図をそれぞれ示す。

先ず、金属部品の成形下死点における応力分布（応力テンソル）を計算する（ステップＳ１０１）。
ここでは、例えば市販の成形シミュレーションプログラム（ＥＳＩ社製；商品名ＰＡＭ−ＳＴＡＭＰ）を用いて、成形下死点における応力分布を計算する。
第１実施形態による成形シミュレーションにおける主要な解析条件を、以下の表１に示す。

ここで、成形シミュレーションによって得られた離型後のスプリングバック変位量の分布を図６に示す。ここでは、変位量を表示濃度に対応して等高線で示しており、表示濃度が濃いほど変位量も大きい。図６に示されるように、特にＡ点におけるＺ軸方向の変位が大きく、金属部品に反り及び捩れが発生していることが判る。

続いて、第１の計算部１０１により、ステップＳ１０１で得られた応力テンソルに基づいて、例えば式（３）の演算を行い、各有限要素（例えば４節点シェル要素）又は複数の有限要素ごとに要素等価節点力ベクトルを計算する（ステップＳ１０２）。
続いて、第２の計算部１０２により、第１の計算部１０１によって計算された１つ又は複数の有限要素ごとの要素等価節点力ベクトルを金属薄板の全領域又は特定領域に対して積分し、その領域の全等価節点力ベクトルを計算する（ステップＳ１０３）。

続いて、特定部１０３により、第２の計算部１０２によって計算された全等価節点力ベクトルの成分のうち、曲げモーメント（Ｍｘ，Ｍｙ，Ｍｚ）について、金属薄板の成形下死点における全領域のうちで、曲げモーメント量の大きい部位の位置を、金属薄板におけるスプリングバックの発生原因部位として特定する（ステップＳ１０４）。

そして、表示部１０４により、曲げモーメント量を、例えば金属薄板の形状に対応させて画像として表示する（ステップＳ１０５）。ここで、ステップＳ１０５をステップＳ１０４の後に行うように記載したが、ステップＳ１０５をステップＳ１０４の先に行うようにしても良い。また、ステップＳ１０４を省略する場合もある。
なお、ステップＳ１０５では、曲げモーメント量を画像表示する代わりに、或いはこれと共に、特定部１０３の計算結果を画像表示するようにしても良い。

表示部１０４による曲げモーメント量の画像表示の一例を図７Ａ及び図７Ｂに示す。図７Ａにｘ軸回りの曲げモーメントＭｘを、図７Ｂにｙ軸回りの曲げモーメントＭｙをそれぞれ示す。ここでは、曲げモーメント量を表示濃度に対応して等高線で示しており、表示濃度が濃いほど曲げモーメント量も大きい。図示のように、金属部品の曲げモーメント量は局所化しており、図７Ａ及び図７Ｂの例ではスプリングバックの発生原因部位として２箇所が明示されている。

ここで、第１実施形態の比較例として、特許文献３の手法により、第１実施形態１と同じ金属部品を用いて成形シミュレーションを行った。
ここでは、ステップＳ１０１により得られた応力分布に基づき、金属部品を複数領域に分割し、応力の高い領域を順次開放し（応力を０にする）、スプリングバック解析（金型拘束を開放し、弾性回復による変形量を計算）を行うことにより、大きな変形量が得られる領域を特定する。

比較例による結果（板厚中央における最大主応力）を図８に示す。ここでは、応力を表示濃度に対応して等高線で示しており、表示濃度が濃いほど応力も大きい。図８に示されるように、大きな応力を示す領域が複数箇所且つ広範囲に分布しており、スプリングバックの発生原因部位を正確に特定することは難しい。また、領域分割の仕方に表示結果が依存するため、正確な変形量表示が困難となる。

これに対して第１実施形態では、図７Ａ及び図７Ｂに示したように、図８とは異なり、例えばｘ軸及びｙ軸まわりの曲げモーメントＭｘ、Ｍｙは、金属部品の稜線の一部に、正負一対で局所的に高い領域がある。従ってこの部分の曲げモーメントがスプリングバック量に対して大きく影響を及ぼすことが容易に特定できる。大きく影響を及ぼす領域が特定できれば、この部分の部品形状を僅かに修正することで、効率良くスプリングバックを低減し、プレス部品の寸法精度を向上させることが容易に可能となる。

以上説明したように、第１実施形態によれば、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速で確実に、金属板の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。

次に、本発明の第２実施形態について、図９〜１４を参照しながら詳細に説明する。

まず、第１実施形態と同様に、金属薄板１０に設定した有限要素について、上記（１）〜（４）式を用いて、全等価節点力ベクトル｛ｆ｝を計算する。ここで、全等価節点力ベクトル｛ｆ｝は、剛性行列（マトリックス）［Ｋ］と、スプリングバックを評価するための特定な節点変位（ｕ_ｉ）とを用いて、以下の式（５）に示すように外力ベクトルとして示される。

続いて、適当な拘束条件を付与した後、剛性マトリックス［Ｋ］の逆行列［Ｋ］^−１を通常の方法で求める。部品全体のスプリングバック量（ｕ）は、通常行われるように、次式のように求まる。

ここで、スプリングバックを評価するための特定な節点変位（ｕ_ｉ）は、次式のように求めることができる。

本発明者らは、式（７）における総和記号中の(k^-1 _ijf_j)は、特定位置のスプリングバック変位量（u_i）に対する、外力ベクトル成分（f_j）毎の寄与（変位）を表していることを発見した。即ち、特定位置のスプリングバック変位と方向が同じ（同符号）で絶対値が大きければ、特定位置のスプリングバック変位に対して正の寄与（スプリングバックを促進）が大であることが判る。また、特定位置のスプリングバック変位と方向が逆（異符号）であれば特定位置のスプリングバック変位に対して負の寄与（スプリングバックを抑制）があることが判る。更に、絶対値が小さければ、特定位置のスプリングバック変位に対しては寄与が小さいことが判る。従って、この節点量を、特定位置の変位量に対する、外力ベクトル成分の寄与として表示することにより、大規模な行列計算を繰り返すことなく、スプリングバック発生の原因となる部位を効率良く特定することができる。

（第２実施形態による成形シミュレーション装置）
図９は、第２実施形態による成形シミュレーション装置の概略構成を示すブロック図である。

この成形シミュレーション装置は、有限要素法を用いて金属薄板の成形下死点における成形シミュレーションを行う。この成形シミュレーション装置は、各要素等価節点力ベクトルを算出する第１の計算部２０１と、全等価節点力ベクトルを算出する第２の計算部２０２と、全体剛性行列の逆行列を計算する第３の計算部２０３と、全等価節点力ベクトルを外力ベクトルとして、全等価節点力ベクトルと全体剛性行列の逆行列とを乗算する第４の計算部２０４と、第４の計算部２０４の計算結果を表示する表示部２０５とを有する。ここで、第１〜第４の計算部２０１〜２０４は、例えばコンピュータの中央処理装置（ＣＰＵ）の各機能として実現される。

第１の計算部２０１は、与えられた応力テンソル（式（２）参照）に基づいて、例えば式（３）の演算を行い、各有限要素（例えば４節点シェル要素）ごと又は複数の有限要素ごとに要素等価節点力ベクトルを計算する。
第２の計算部２０２は、第１の計算部２０１により計算された各有限要素ごと（又は複数の有限要素ごと）の要素等価節点力ベクトルを金属薄板の全領域に対して積分し（式（４）参照）、全等価節点力ベクトルを計算する。

第３の計算部２０３は、全体剛性マトリクスの逆行列を計算する。全体剛性マトリクスは、例えば非特許文献１に示されるように、前述した各種有限要素に対応する変位−歪み関係マトリックス、一般的な線形弾性構成則に基づく応力−歪み関係マトリックス、及び座標変換マトリックスを用いて求められる。逆行列の計算法としては、通常行われる手法を用いれば良い。

第４の計算部２０４は、全等価節点力ベクトルを外力ベクトルとして、全体剛性マトリクスの逆行列と全等価節点力ベクトルとを乗算する（式（５）〜（７）参照）。この乗算では、式（７）における総和記号中の(k^-1 _ijf_j)が外力ベクトルの成分ごとの寄与度を示している。この乗算により得られた特定位置のスプリングバック変位量は、各寄与度の総和を表している。

表示部２０５は、第４の計算部２０４により得られた、金属薄板の特定位置の変位量に対する外力ベクトルの成分ごとの寄与度を表示する。具体的には、例えば金属薄板の特定位置に対して、外力ベクトルの成分ごとの寄与度に応じて色付けし、当該特定位置に対するスプリングバックの発生原因部位を視認できるように表示する。

（第２実施形態による成形シミュレーション方法）
以下、上記の成形シミュレーション装置を用いた成形シミュレーション方法について説明する。

図１０は、第２実施形態による成形シミュレーション方法をステップ順に示すフロー図である。
金属薄板としては、ここでは図１１Ａ及び図１１Ｂに示す金属部品（高張力鋼板）を成形シミュレーションの対象とする。ここで、図１１Ａが金属部品の全体図、図１１Ｂが部分拡大図をそれぞれ示す。

先ず、金属部品の成形下死点における応力分布（応力テンソル）を計算する（ステップＳ２０１）。
ここでは、例えば市販の成形シミュレーションプログラム（ＥＳＩ社製；商品名ＰＡＭ−ＳＴＡＭＰ）を用いて、成形下死点における応力分布を計算する。
第２実施形態による成形シミュレーションにおける主要な解析条件は、第１実施形態と同じく、表１で示される条件を用いた。

ここで、成形シミュレーションによって得られた離型後のスプリングバック変位量の分布を図１２に示す。ここでは、変位量を表示濃度に対応して等高線で示しており、表示濃度が濃いほど変位量も大きい。図１２に示されるように、特にＡ点におけるＺ軸方向の変位が大きく、金属部品に反り及び捩れが発生していることが判る。

続いて、第１の計算部２０１により、ステップＳ２０１で得られた応力テンソルに基づいて、例えば式（３）の演算を行い、各有限要素（例えば４節点シェル要素）ごと又は複数の有限要素ごとに要素等価節点力ベクトルを計算する（ステップＳ２０２）。
続いて、第２の計算部２０２により、第１の計算部２０１によって計算された各有限要素ごと（又は複数の有限要素ごと）の要素等価節点力ベクトルを金属薄板の全領域（又は特定の領域）に対して積分し、全等価節点力ベクトルを計算する（ステップＳ２０３）。

続いて、第３の計算部２０３により、全体剛性マトリクスの逆行列を計算する（ステップＳ２０４）。

続いて、第４の計算部２０４により、全等価節点力ベクトルを外力ベクトルとして、全体剛性マトリクスの逆行列と全等価節点力ベクトルとを乗算する（ステップＳ２０５）。
そして、表示部２０５により、金属薄板の所定（特定）位置に対するスプリングバック変位量の外力ベクトルの成分ごとの寄与度を画像表示する（ステップＳ２０６）。

表示部２０５による、金属薄板の特定位置に対するスプリングバック変位量の外力ベクトルの成分ごとの寄与度の一例を図１３Ａ及び図１３Ｂに示す。ここでは、図１２のＡ点のＺ軸方向変位に対する外力ベクトルの成分である曲げモーメント量（Ｍｘ，Ｍｙ）の寄与度（Ｄｚ）分布を等高線で示す。ここで、寄与（変位）は、Ａ点のｚ方向変位で除することで無次元化してある。即ち図示の例では、Ａ点における変位に対する外力ベクトル成分の寄与率を表している。図１３Ａにｘ軸回りの曲げモーメントＭｘのＤｚ分布を示す。図１３Ｂにｙ軸回りの曲げモーメントＭｙのＤｚ分布を示す。図１３Ａ及び図１３Ｂに示される例では、Ｄｚを表示濃度に対応して等高線で示しており、表示濃度が濃いほどＤｚも大きい。

ここで、第２実施形態の比較例として、特許文献３の手法により、第２実施形態と同じ金属部品を用いて成形シミュレーションを行った。
ここでは、ステップＳ２０１により得られた応力分布に基づき、金属部品を複数領域に分割し、応力の高い領域を順次開放し（応力を０にする）、スプリングバック解析（金型拘束を開放し、弾性回復による変形量を計算）を行うことにより、大きな変形量が得られる領域を特定する。

比較例による結果（板厚中央における最大主応力）を図１４に示す。ここでは、応力を表示濃度に対応して等高線で示しており、表示濃度が濃いほど応力も大きい。図１４に示されるように、大きな応力を示す領域は複数箇所且つ広範囲に分布している。従って、スプリングバックの発生原因部位を正確に特定することは難しい。また、領域分割の仕方に表示結果が依存するため、正確な変形量表示が困難となる。

これに対して第２実施形態では、図１３Ａ及び図１３Ｂに示したように、図１４とは異なり、図１２のＡ点におけるＤｚ分布の影響の大きさを、ピンポイント且つ定量的に把握することができる。影響の大きな領域が特定できれば、この部分の部品形状を僅かに修正することで、効率良くスプリングバックを低減し、プレス部品の寸法精度を向上させることが容易に可能となる。

以上説明したように、第２実施形態によれば、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速で確実に、金属薄板の成形時におけるスプリングバックの発生部位を効率良く特定し、正確な成形に供することが可能となる。

次に、本発明の第３実施形態について、図１５〜２０を参照しながら詳細に説明する。

まず、第２実施形態と同様に、金属薄板１０に設定した有限要素について、上記（１）〜（７）式を用いて、スプリングバックを評価するための特定の節点変位（ｕ_ｉ）を計算する。ここで、計算されたスプリングバックを評価するための特定な節点変位を、基準となる第一の変位（ｕ_i ⁽¹⁾）とする。

続いて、各要素毎に計算された要素等価節点力ベクトル｛ｆ｝_ｅを、全等価節点力ベクトル｛ｆ｝から差し引いた修正全等価節点力ベクトル｛ｆ’｝を計算する。

ここで、要素等価節点力ベクトル｛ｆ｝_ｅは、前述した（３）式を用いて計算しても良いし、あるいは特定な成分、例えば面内力やモーメントのみを用いても良い。また、後述するようにスプリングバックに対する要素の寄与度を求めることが目的であるため、ある係数でスカラー倍した値を用いても良い。

ついで、式（６）で求めた逆行列［Ｋ］^−１と修正全等価節点力ベクトル｛ｆ’｝を用いて、スプリングバックを評価するための特定の節点に対する第二の変位（ｕ_i ⁽²⁾）を式（９）のように求める。

式（８）と式（９）の計算を、要素毎に行い、第一の変位（ｕ_i ⁽¹⁾）と第二の変位（ｕ_i ⁽²⁾）の変化量を、例えば以下のように求める。

式（１１）に示すように、第一の変位（ｕ_i ^(１)）の絶対値よりも第二の変位（ｕ_i ⁽²⁾）の絶対値が小さければ、当該要素の残留応力は、開放することでスプリングバックを低減させる作用があることを示す。
また、式（１２）に示すように、第一の変位（ｕ_i ^(１)）の絶対値よりも第二の変位（ｕ_i ⁽²⁾）の絶対値が大きければ、当該要素の残留応力は、スプリングバックを増大させる作用があることを示す。
更に、式（１３）に示すように、第一の変位（ｕ_i ^(１)）の絶対値と第二の変位（ｕ_i ⁽²⁾）の絶対値とがほぼ等しければ、当該要素の残留応力は、スプリングバック寄与していないことを示す。

また、式（１０）で計算される第一の変位（ｕ_i ^(１)）と第二の変位（ｕ_i ⁽²⁾）との変化量を要素ごとに表示することにより、スプリングバック発生の原因となる部位を効率良く特定することができる。

さらに、式（８）および式（９）の計算は、大規模な行列反転操作を含まない単純なベクトル計算である。従って計算負荷は小さく、前述した要素毎の変位量計算を全要素について行っても、短時間で計算可能である。

ここで、第一の変位（ｕ_i ^(１)）と第二の変位（ｕ_i ⁽²⁾）の変化量の算出方式は、式（１０）に限定されるものでない。これに替わって、例えば第一の変位で無次元化した、下記式（１４）等の様な変化量の算出方式が考えられる。

さらに例えば、特定の節点変位は、単一の節点でなく、下記式（１５）等の様な複数の節点変位量の算術式で表されても良い。

（第３実施形態による成形シミュレーション装置）
図１５は、第３実施形態による成形シミュレーション装置の概略構成を示すブロック図である。

この成形シミュレーション装置は、有限要素法を用いて金属薄板の成形下死点における成形シミュレーションを行うものである。この成形シミュレーション装置は、各要素等価節点力ベクトル｛ｆ｝_ｅを算出する第１の計算部３０１と、全等価節点力ベクトル｛ｆ｝を算出する第２の計算部３０２と、全体剛性行列の逆行列[Ｋ]^−１を計算する第３の計算部３０３と、全等価節点力ベクトル｛ｆ｝を外力ベクトルとして、全等価節点力ベクトル｛ｆ｝と全体剛性行列の逆行列[Ｋ]^−１とを乗算する第４の計算部３０４と、一つの又は複数の要素毎に、修正全等価節点力ベクトル｛ｆ’｝と全体剛性行列の逆行列[Ｋ]^−１とを乗算する第５の計算部３０５と、一つの又は複数の要素毎に、計算部３０４の計算結果と計算部３０５の計算結果との変化量を計算する第６の計算部３０６と、計算部３０６の計算結果を表示する表示部３０７とを有する。ここで、第３０１〜第３０６の計算部３０１〜３０６は、例えばコンピュータの中央処理装置（ＣＰＵ）の各機能として実現される。

第１の計算部３０１は、与えられた応力テンソル（式（２）参照）に基づいて、例えば式（３）の演算を行い、各有限要素（例えば４節点シェル要素）ごと又は複数の有限要素ごとに要素等価節点力ベクトル｛ｆ｝_ｅを計算する。

第２の計算部３０２は、第１の計算部３０１により計算された各有限要素ごと（又は複数の有限要素ごと）の要素等価節点力ベクトル｛ｆ｝_ｅを金属薄板の全領域（又は特定領域）に対して積分し（式（４）参照）、その領域の全等価節点力ベクトル｛ｆ｝を計算する。

第３の計算部３０３は、全体剛性マトリクスの逆行列[Ｋ]^−１を計算する。全体剛性マトリクスは、例えば非特許文献１に示されるように、前述した各種有限要素に対応する変位−歪み関係マトリックス、一般的な線形弾性構成則に基づく応力−歪み関係マトリックス、及び座標変換マトリックスを用いて求められる。逆行列[Ｋ]^−１の計算法としては、通常行われる手法を用いれば良い。

第４の計算部３０４は、全等価節点力ベクトル｛ｆ｝を外力ベクトルとして、全体剛性マトリクスの逆行列[Ｋ]^−１と全等価節点力ベクトル｛ｆ｝とを乗算する（式（５）〜（７）参照）。

第５の計算部３０５は、各要素毎に計算された要素等価節点力ベクトル｛ｆ｝_ｅを、全等価節点力ベクトル｛ｆ｝から差し引いた修正全等価節点力ベクトル｛ｆ’｝を外力として、全体剛性マトリクスの逆行列[Ｋ]^−１と修正全等価節点力ベクトル｛ｆ’｝とを乗算する（式（８）〜（９）参照）。

第６の計算部３０６は、第４の計算部３０４から求められた第一の変位（ｕ_i ⁽¹⁾）と、第５の計算部３０５から求められた第二の変位（ｕ_i ^(２)）を用いて、スプリングバックを評価するための特定の節点に対する変位の変化量を求める。（式（１０）参照）。

表示部３０７は、第６の計算部３０６により得られた、金属薄板の特定位置の変位量に対する要素ごとの寄与度を表示するものである。具体的には、例えば金属薄板の特定位置に対して、要素ごとの寄与度に応じて色付けし、当該特定位置に対するスプリングバックの発生原因部位を視認できるように表示する。

（第３の実施形態による成形シミュレーション方法）
以下、上記の成形シミュレーション装置を用いた成形シミュレーション方法について説明する。

図１６は、第３実施形態による成形シミュレーション方法をステップ順に示すフロー図である。金属薄板としては、ここでは図１７Ａ及び図１７Ｂに示す金属部品（高張力鋼板）を成形シミュレーションの対象とする。ここで、図１７Ａが金属部品の全体図、図１７Ｂが部分拡大図をそれぞれ示す。

先ず、金属部品の成形下死点における応力分布（応力テンソル）を計算する（ステップＳ３０１）。
ここでは、例えば市販の成形シミュレーションプログラム（ＥＳＩ社製；商品名ＰＡＭ−ＳＴＡＭＰ）を用いて、成形下死点における応力分布を計算する。
第３実施形態による成形シミュレーションにおける主要な解析条件は、第１実施形態と同じく、表１で示される条件を用いた。

ここで、成形シミュレーションによって得られた離型後のスプリングバック変位量の分布を図１８に示す。ここでは、変位量を表示濃度に対応して等高線で示しており、表示濃度が濃いほど変位量も大きい。図１８に示されるように、特にＡ点におけるＺ軸方向の変位が大きく、金属部品に反り及び捩れが発生していることが判る。

続いて、第１の計算部３０１により、ステップＳ３０１で得られた応力テンソルに基づいて、例えば式（３）の演算を行い、各有限要素（例えば４節点シェル要素）ごと又は複数の有限要素ごとに要素等価節点力ベクトル｛ｆ｝_ｅを計算する（ステップＳ３０２）。

続いて、第２の計算部３０２により、第１の計算部３０１によって計算された各有限要素ごと（又は複数の有限要素ごと）の要素等価節点力ベクトル｛ｆ｝_ｅを金属薄板の全領域（又は特定の領域）に対して積分し、その領域の全等価節点力ベクトル｛ｆ｝を計算する（ステップＳ３０３）。

続いて、第３の計算部３０３により、全体剛性マトリクスの逆行列［Ｋ］^−１を計算する（ステップＳ３０４）。

続いて、第４の計算部３０４により、全等価節点力ベクトル｛ｆ｝を第１の外力ベクトルとして、全体剛性マトリクスの逆行列［Ｋ］^−１と全等価節点力ベクトル｛ｆ｝とを乗算する（ステップＳ３０５）。

続いて、第５の計算部３０５は、各要素毎に計算された要素等価節点力ベクトル｛ｆ｝_ｅを、全等価節点力ベクトル｛ｆ｝から差し引いた修正全等価節点力ベクトル｛ｆ’｝を第２の外力として、全体剛性マトリクスの逆行列［Ｋ］^−１と修正全等価節点力ベクトル｛ｆ’｝とを乗算する（ステップＳ３０６）。

続いて、第６の計算部３０６は、第４の計算部３０４から求められた第一の変位（ｕ_i ⁽¹⁾）と、第５の計算部３０５から求められた第二の変位（ｕ_i ^(２)）とを用いて、スプリングバックを評価するための特定な節点に対する変位の変化量、すなわち要素ごとの寄与度を求める。（ステップＳ３０７）。

そして、表示部３０７により、金属薄板の所定（特定）位置に対するスプリングバック変位量の要素ごとの寄与度を画像表示する（ステップＳ３０８）。

表示部３０７による、金属薄板の特定位置に対するスプリングバック変位量の要素ごとの寄与度の一例を図１９に示す。図１９における一点鎖線は、金属薄板の成形下死点における形状を示す輪郭線である。図１９では、図１８のＡ点におけるｚ方向変位に対する、要素の寄与度分布を等高線で示している。ここで、寄与度は、Ａ点のｚ方向変位で除することで無次元化しており、表示濃度が濃いほど寄与度も大きい。

ここで、第３実施形態の比較例として、特許文献１の手法により、第３実施形態と同じ金属部品を用いて成形シミュレーションを行った。
ここでは、ステップＳ３０１により得られた応力分布に基づき、金属部品を複数領域に分割し、応力の高い領域を順次開放し（応力を０にする）、スプリングバック解析（金型拘束を開放し、弾性回復による変形量を計算）を行うことにより、大きな変形量が得られる領域を特定する。

比較例による結果（領域ごとの寄与度）を図２０に示す。ここでは、寄与度を表示濃度に対応して等高線で示しており、表示濃度が濃いほど寄与度も大きい。図２０に示されるように、寄与度の大きな領域は複数箇所且つ広範囲に分布しており、スプリングバックの発生原因部位を正確に特定することは難しい。また、領域分割の仕方に表示結果が依存するため、正確な変形量表示が困難となる。

これに対して第３実施形態では、図２０とは異なり、図１９に示すように、図１８のＡ点変位に対する影響の大きさを、ピンポイント且つ定量的に把握することができる。影響の大きな領域がピンポイントに特定できれば、この部分の部品形状を僅かに修正することで、効率良くスプリングバックを低減し、プレス部品の寸法精度を向上させることが容易に可能となる。

以上説明したように、第３実施形態によれば、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速で確実に、金属薄板の成形時におけるスプリングバックの発生部位を効率良く特定し、正確な成形に供することが可能となる。

以上説明してきた、各実施形態による計算時間の比較を表２に示す。ここで計算時間は比較例を１００とした場合の相対値である。

次に、本発明の第４実施形態について、図２１を参照しながら詳細に説明する。

上述した第１〜３実施形態による成形シミュレーション装置を構成する各構成要素（表示部１０４、２０５、３０７を除く）等の機能は、コンピュータのＲＡＭやＲＯＭ等に記憶されたプログラムが動作することによって実現できる。同様に、成形シミュレーション方法の各ステップ（図４のステップＳ１０１〜Ｓ１０５、図１０のＳ２０１〜Ｓ２０６、図１６のステップＳ３０１〜Ｓ３０８等）は、コンピュータのＲＡＭやＲＯＭ等に記憶されたプログラムが動作することによって実現できる。このプログラム及び当該プログラムを記録したコンピュータ読み取り可能な記憶媒体は本発明に含まれる。

具体的に、前記プログラムは、例えばＣＤ−ＲＯＭのような記録媒体に記録し、或いは各種伝送媒体を介し、コンピュータに提供される。前記プログラムを記録する記録媒体としては、ＣＤ−ＲＯＭ以外に、フレキシブルディスク、ハードディスク、磁気テープ、光磁気ディスク、不揮発性メモリカード等を用いることができる。他方、前記プログラムの伝送媒体としては、プログラム情報を搬送波として伝搬させて供給するためのコンピュータネットワークシステムにおける通信媒体を用いることができる。ここで、コンピュータネットワークとは、ＬＡＮ、インターネットの等のＷＡＮ、無線通信ネットワーク等であり、通信媒体とは、光ファイバ等の有線回線や無線回線等である。

また、本発明に含まれるプログラムとしては、供給されたプログラムをコンピュータが実行することにより上述の実施形態の機能が実現されるようなもののみではない。例えば、そのプログラムがコンピュータにおいて稼働しているＯＳ（オペレーティングシステム）或いは他のアプリケーションソフト等と共同して上述の実施形態の機能が実現される場合にも、かかるプログラムは本発明に含まれる。また、供給されたプログラムの処理の全て或いは一部がコンピュータの機能拡張ボードや機能拡張ユニットにより行われて上述の実施形態の機能が実現される場合にも、かかるプログラムは本発明に含まれる。

例えば、図２１は、パーソナルユーザ端末装置の内部構成を示す模式図である。この図２１において、４００はＣＰＵ４０１を備えたパーソナルコンピュータ（ＰＣ）である。ＰＣ４００は、ＲＯＭ４０２またはハードディスク（ＨＤ）４１１に記憶された、又はフレキシブルディスクドライブ（ＦＤ）４１２より供給されるデバイス制御ソフトウェアを実行する。このＰＣ４００は、システムバス４０４に接続される各デバイスを総括的に制御する。

ＰＣ４００のＣＰＵ４０１、ＲＯＭ４０２またはハードディスク（ＨＤ）４１１に記憶されたプログラムにより、第１実施形態の図４におけるステップＳ１０１〜Ｓ１０５、第２実施形態の図１０におけるステップＳ２０１〜Ｓ２０６、及び第３実施形態の図１６におけるステップＳ３０１〜Ｓ３０８の手順等が実現される。

４０３はＲＡＭであり、ＣＰＵ４０１の主メモリ、ワークエリア等として機能する。４０５はキーボードコントローラ（ＫＢＣ）であり、キーボード（ＫＢ）４０９や不図示のデバイス等からの指示入力を制御する。

４０６はＣＲＴコントローラ（ＣＲＴＣ）であり、ＣＲＴディスプレイ（ＣＲＴ）４１０の表示を制御する。４０７はディスクコントローラ（ＤＫＣ）である。ＤＫＣ４０７は、ブートプログラム、複数のアプリケーション、編集ファイル、ユーザファイルそしてネットワーク管理プログラム等を記憶するハードディスク（ＨＤ）４１１、及びフレキシブルディスク（ＦＤ）４１２とのアクセスを制御する。ここで、ブートプログラムとは、起動プログラム：パソコンのハードやソフトの実行（動作）を開始するプログラムである。

４０８はネットワーク・インターフェースカード（ＮＩＣ）で、ＬＡＮ４２０を介して、ネットワークプリンタ、他のネットワーク機器、或は他のＰＣと双方向のデータのやり取りを行う。

上記のパーソナルユーザ端末装置によれば、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速で確実に、金属板の成形時におけるスプリングバックの発生部位を効率良く特定し、正確な成形に供することができる。

次に、本発明の第５実施形態について、図２２〜３０を参照しながら詳細に説明する。

図２２は、本発明のプレス成形時におけるスプリングバック抑制の対象となる金属薄板のプレス部品（以下本部品という）の形状の一例を示す概略斜視図である。
本部品５０１を、薄鋼板、アルミ薄板などの金属薄板、あるいはプラスチックや複合材料等の弾塑性材料を素材として、上型および下型からなる金型を用いて、冷間プレス成形すると、離型後に、弾性回復（スプリングバック）により、部品形状は、成形下死点での形状（即ち成形目標形状）から変化する。そのため、所定の部品寸法精度が得られず、溶接不良、組立精度悪化などの不具合が発生する。

図２３に、本部品５０１のスプリングバックの一例をスプリングバック変位量の分布コンタ図で示す。色の濃淡は、成形下死点からのZ軸（紙面直角方向）変位量を示す。図２３に示す本部品５０１のＡ、Ｂ、Ｃ、Ｄ点をつないでなる平面は、Ｘ−Ｙ平面（紙面水平方向）に平行に成形されるべきであるが、湾曲形状の外側端点（Ａ、Ｂ点）がＺ軸方向に持ち上がり断面が捩れていることを示している。
前記スプリングバックを合理的に抑制することが、本発明の目的である。

スプリングバックを抑制するには、部品自体の弾性剛性を向上させる（第１の対策方法）、あるいはスプリングバックを誘起する内部残留応力を低減させる（第２の対策方法）、ということが一般的に行われる。第１の対策方法に対しては、部品形状の変更や形状凍結ビードの付与などが一般に用いられる。また第２の対策方法に対しては、ビード付与、多段成形や温間、熱間成形などがある。
これらの対策によっても、寸法精度が不十分であれば、スプリングバック量を予め見込んで、金型形状を製品と異なる形状に修正することが通常行われる（第３の対策方法）。しかしながら、見込み量の精度は様々な因子に敏感であり、金型修正コストも増大する。

一方、部品形状が複雑になると、成形下死点での内部残留応力は一様でなく複雑に分布する。図２４に、本部品５０１の成形下死点での内部残留応力を、有限要素法に基づくプレス成形シミュレーションによって計算予測した一例を示す。図２４においては板厚中心での主応力分布を色の濃淡で表しており、色の濃い部分で示される内部残留応力の大きな部位が複数あり、また引張、圧縮が混在している。
プレス成形シミュレーションとしては上記第１実施形態から第４実施形態のシミュレーションをもちいることができる。

本発明のスプリングバック抑制方法では、本部品１のスプリングバックを誘起する部位を特定する第一の工程と、スプリングバックを誘起する部位に予め、1個以上のエンボスを成形する第二の工程と、所定の部品形状になるようにエンボスを平坦に潰す第三の工程とを行う。
本発明の、第一の工程では、複雑な内部残留応力分布を示す複雑形状の本部品１に対して、本部品１の内部応力を部分的に変化させることでスプリングバックが変化する領域を、スプリングバックを誘起している領域（部位）として抽出し、特定する。

抽出方法として、有限要素法に基づくプレス成形シミュレーションによって計算予測した本部品１の内部残留応力を、部分領域毎に変化させて再度計算予測し、内部残留応力の変化の前後でスプリングバック量が大きく変わる領域を、スプリングバックを誘起する部位として特定する。内部残留応力の変化のさせ方は、部分領域の応力成分を板厚方向全積分点（応力評価点）で一律０（開放）（成形下死点の応力を０とする）しても良い。あるいは、膜応力成分（板厚中心積分点の応力）と曲げ応力成分（板厚方向積分点毎応力から膜応力を差し引いた）に分解し、それぞれを０とすることで、面内応力と曲げ応力の影響を分離して評価することができる。また、部分領域の大きさは、フランジ部、曲げＲ部、ウェブ部など形状の特徴に合わせた領域分割でもよい。しかし、望ましくは、前記形状を表現するに十分細かい有限要素メッシュ分割をそのまま部分領域とすることで、高い分解能で抽出が可能となる。

図２５は、図２４に示す有限要素法に基づくプレス成形シミュレーションにおける有限要素に対応する部分領域ごとに内部残留応力を開放し、開放前後のスプリングバック量の差を、開放前のスプリングバック量で基準化した値の分布を色の濃淡で表したコンタ図である。図２５において濃い色で示されるように、本部品５０１の一部である変化の大きい部位（開放することでスプリングバックが低減した領域）Ｅ、Ｆを、スプリングバックを誘起する部位Ｅ、Ｆとして特定できる。図２５に示すスプリングバックを誘起する部位Ｅ、Ｆは、図２４に示した内部残留応力の大きな部位とは異なり、必ずしも応力の大きい部位がスプリングバックを誘起しているわけでは無いことを示している。

また、本部品５０１の内部応力を部分的に変化させることでスプリングバックが変化する領域を、スプリングバックを誘起する部位として特定する別の抽出方法もある。すなわち、実物のプレス部品を用いて、一部を切り欠いて、あるいは穴を開けるなどして、部分的に応力を開放し、その前後でスプリングバック量が大きく変化する部位を、スプリングバックを誘起する部位として特定してもよい。

続いて、本発明の第二の工程では、前記方法にて特定されたスプリングバックを誘起する部位に、プレス第一工程でエンボスを形成する。図２６には、プレス第一工程において成形された本部品５０１の平面図を示す。ここでは、図２５に示したスプリングバックを誘起する部位Ｅ、Ｆに、同一寸法の円形エンボス５０２、５０３を２個成形している。
図２６におけるエンボス部Ａ−Ａ断面の成形された本部品５０１および金型（下型５０４、上型５０５）配置を図２７に示す。図２７中、エンボス成形部は、円筒状工具５０６を、下型５０４に挿入する構造としている。工具５０６として下型５０４の凹部５０４ａに勘合する高さ寸法の異なる複数の着脱式の円筒状工具５０６を用い、円筒状工具５０６を交換することで、容易にエンボス５０２、５０３の高さ調整が行える。エンボス形状（工具５０６の平面形状）は、円形に限定されるものでなく、楕円、四角形など任意の形状が可能だが、円筒形状とすることで、加工が容易で、種々の高さの工具５０６をエンボス５０２、５０３の高さ寸法調整用に準備しても、低コスト化が図れる。

続いて、本発明の第三の工程、すなわちプレス第二工程では、前記プレス第一工程でエンボスを成形された成形部品を、所定の部品形状になるようにエンボスの無い平坦な工具で押しつぶす。図２８には、プレス第二工程において成形された本部品５０１の平面図を示す。また図２７に対応するプレス第二工程Ａ−Ａ断面の成形された本部品５０１および金型（下型５０７、上型５０８）配置を図２９に示す。一旦張出したエンボスを平坦に押し潰すことにより、繰り返し曲げ−曲げ戻しの効果で、スプリングバックを誘起する内部残留応力を緩和することができ、スプリングバック量が大幅に低減する。

また、金型の調整段階では、プレス第一工程で、前記種々の高さの円筒工具でエンボス高さを調整する。続いてプレス第二工程でエンボスを潰し平坦にする。これにより、製品形状を変更することなく、スプリングバックを容易に制御、抑制することができる。また、図２７に示す前記円筒工具６の端部を湾曲したエンボス形状でなく平坦にすることで、不要なエンボスを削除することも可能である。また、同じく図２７に示す下型５０４のエンボス凸部として凹部５０４ａに前記端部平坦な円筒工具を挿入し、さらに上型５０５のエンボス凹部５０５ａを着脱式の工具で平坦に埋めることで、プレス第二工程を、別の金型を使わずに、プレス第一工程と同一金型で成形することも可能である。

図２７および図２８では、プレス第一工程とプレス第二工程で、エンボス部以外は同一形状としているが、これに限定されるものでなく、角部曲率半径や他部品取り付け座面など種々の形状を、工程ごとに成形しても構わない。

図２２は、本発明を適用した、プレス部品（本部品５０１）の斜視図である。本部品５０１の材料は、板厚１．２ｍｍ、引張強度９８０ＭＰａ級の冷延高張力鋼板である。プレス形成用金型を製作するに先立ち、成形シミュレーションプログラム（ＥＳＩ社製；商品名ＰＡＭ−ＳＴＡＭＰ）を用いて、成形下死点における応力分布に基づき、スプリングバックを誘起している領域を抽出した。
第５実施形態による成形シミュレーションにおける主要な解析条件は、第１実施形態と同じく、表１で示される条件を用いた。

ここで、成形シミュレーションによって得られた離型後の本部品１のスプリングバック変位量の分布を図２３に示す。ここでは、変位量を表示濃度に対応して等高線で示しており、表示濃度が濃いほど変位量も大きい。図示のように、特にＡ点およびＢ点におけるＺ軸方向（紙面垂直方向）の変位が大きく、本部品５０１に反り及び捩れが発生していることが判る。また、図２４には、成形下死点における最大主応力分布コンタ図を示す。

続いて、スプリングバックを誘起している領域を特定するため、部分領域毎に、成形下死点の応力を強制的に０として、離型後のスプリングバック変位量を計算した。ここで、捩れに対するスプリングバック量を表す指標として、Ａ点およびＣ点のＺ軸方向変位ＤＡ、ＤＣの差、ＤＡ−ＤＣ（＝Ｄ_Ａ−Ｃ）を用い、応力開放前後のスプリングバック量から、その領域のスプリングバック量に対する影響度を、式（１６）で表した。

影響度＝Ｄ_Ａ−Ｃ（応力開放後）／Ｄ_Ａ−Ｃ（応力開放前）−１ ・・・（１６）

部分領域を有限要素単位に一致させ、要素ごとに（１６）式の計算を行った結果、影響度の分布コンタ図を図２５に示す。影響度の分布を色の濃淡で表しており、本部品５０１の一部に変化の大きい部位Ｅ、Ｆ、すなわちスプリングバックを誘起する部位Ｅ、Ｆが特定できる。また、スプリングバックを誘起する部位Ｅ、Ｆは、図２４に示した内部残留応力の分布とは異なり、必ずしも応力の大きい部位がスプリングバックを誘起しているのでは無いことが明らかとなった。

続いて、前記スプリングバックを誘起している領域Ｅ、Ｆを成形シミュレーションによって特定した結果に基づき、プレス成形用金型を試作し、成形実験を行った。
前記影響度の大きいウェブ面の部位Ｅ、Ｆの２箇所（Ｅ点、Ｆ点）に、プレス第一工程でエンボス５０２、５０３を形成した。図２６にプレス第一工程の成形部品の平面図を示す。

エンボス５０２、５０３はいずれも直径φＤ＝１０ｍｍの球殻形状で、成形用の金型は、図３０に示すような高さＨ＝０、１、２、３ｍｍの４種類の円筒工具５０６を準備し、図２７に示すように、下型５０４に挿入する構造とした。
前記２箇所Ｅ点およびＦ点のエンボス高さを種々変えたプレス第一工程サンプルを試作し、続いて、図２９に示すような平滑金型（下型５０７、上型５０８）でエンボス５０２、５０３を潰すプレス第二工程で、図２８に示すような製品サンプル（本部品１）を得た。
各製品サンプルに対して、前記スプリングバック量Ｄ_Ａ−Ｃを測定した結果を表３に示す。

エンボス高さ０の比較例に対して、エンボス高さを高くするとスプリングバック量（Ｄ_Ａ−Ｃ）は低減すること、さらに高さを種々調節することで、スプリングバック量をスプリングゴー（Ｄ_Ａ−Ｃが−である（スプリングイン））まで含めて、可変制御可能であることが確認された。
尚、スプリングバックを誘起する部位を特定するために、本発明の成形シミュレーション方法を利用することが好ましい。

以上説明したように、第５実施形態によれば、複雑形状のプレス部品であっても、プレス成形時におけるスプリングバックを誘起する部分をピンポイントに特定し、製品形状を変更することなく、調整容易な金型を用いて、寸法精度に優れたプレス成形品を効率よく得ることが可能となる。以上、金属薄板のスプリングバックを抑制する方法に関して本実施形態を説明してきた。しかし本発明は金属材料に限られるものではなく、プラスチック材や複合材料等の弾塑性材料を成形する場合も含む。また、薄板に限られるものでもなく、線材や、ある程度厚みのある材料を成形する場合も含む。

本発明によれば、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速で確実に、弾塑性材料の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。従って、成形シミュレーションに要する時間の短縮化、及び、成形品の高精度化に寄与することが出来る。

１０１ 第１の計算部
１０２ 第２の計算部
１０３ 特定部
１０４ 表示部
２０１ 第１の計算部
２０２ 第２の計算部
２０３ 第３の計算部
２０４ 第４の計算部
２０５ 表示部
３０１ 第１の計算部
３０２ 第２の計算部
３０３ 第３の計算部
３０４ 第４の計算部
３０５ 第５の計算部
３０６ 第６の計算部
３０７ 表示部
５０１ プレス部品
５０２ エンボス（Ｅ点）
５０３ エンボス（Ｆ点）
５０４ プレス第一工程下型
５０５ プレス第一工程上型
５０６ エンボス成形用工具
５０７ プレス第二工程下型
５０８ プレス第二工程上型

上記の目的を達成するために、本発明は以下の手段を採用した。
（１）本発明の弾塑性材料の成形シミュレーション方法は、第１計算部が、有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する工程と；第２計算部が、計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する工程と；特定部が、計算された前記全等価節点力ベクトルに基づいて、前記弾塑性材料の全領域又は特定の領域のうちで、前記全等価節点力ベクトルの構成成分である曲げモーメントが大きい部位をスプリングバックの発生原因部位として特定する工程と；を含む。
（２）前記（１）に記載の成形シミュレーション方法は、第３の計算部が、全体剛性行列の逆行列を計算する工程と；第４の計算部が、前記全等価節点力ベクトルを外力ベクトルとして、前記全体剛性行列の逆行列と前記全等価節点力ベクトルとを乗算し、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの寄与度を計算する工程と；を更に含んでも良い。
（３）前記（２）に記載のシミュレーション方法は、表示部が、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの前記寄与度を表示する工程を更に含んでも良い。
（４）前記（１）に記載の成形シミュレーション方法は、第３計算部が、全体剛性行列の逆行列を計算する工程と；第４計算部が、前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する工程と；第５計算部が、前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する工程と；第６計算部が、前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する工程と；を更に含んでも良い。
（５）前記（４）に記載の成形シミュレーション方法は、表示部が、前記第１の変位量と前記第２の変位量との間の変化量を表示する工程を更に含んでも良い。
（６）本発明の弾塑性材料の成形方法は、前記（１）に記載のシミュレーション方法により、前記弾塑性材料のスプリングバック発生原因部位を特定する工程と；、前記弾塑性材料の前記スプリングバック発生原因部位に予め、エンボスを成形する工程と；前記エンボスに対し、圧縮応力を与えるように塑性変形させる工程と；を含む。
（７）前記（６）に記載の弾塑性材料の成形方法は、前記エンボス成形工程で成形する前記エンボスが、２個以上の略同一寸法の円形エンボスであっても良い。
（８）前記（６）に記載の弾塑性材料の成形方法は、前記エンボス成形工程が、金型に設けられた凹部に勘合する高さ寸法の異なる複数の着脱式の工具を用いてエンボス形状の高さ寸法を調整する工程を有しても良い。
（９）本発明のプログラムは、コンピュータにより、弾塑性材料の成形シミュレーションを行うに際して、有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する工程と；計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する工程と；計算された前記全等価節点力ベクトルに基づいて、前記弾塑性材料の全領域又は特定の領域のうちで、前記全等価節点力ベクトルの構成成分である曲げモーメントが大きい部位をの大きい部位をスプリングバックの発生原因部位として特定する工程と；をコンピュータに実行させる。
（１０）前記（９）に記載のプログラムは、全体剛性行列の逆行列を計算する工程と；前記全等価節点力ベクトルを外力ベクトルとして、前記全体剛性行列の逆行列と前記全等価節点力ベクトルとを乗算し、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの寄与度を計算する工程と；を更にコンピュータに実行させてもよい。
（１１）前記（９）に記載のプログラムは、全体剛性行列の逆行列を計算する工程と；前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する工程と；前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する工程と；前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する工程と；を更にコンピュータに実行させても良い。
（１２）本発明のコンピュータで読み取り可能な記録媒体は、前記（９）に記載のプログラムを記録している。
（１３）本発明の、弾塑性材料の成形シミュレーションを行う成形シミュレーション装置は、有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する第１の計算部と；計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する第２の計算部と；全体剛性行列の逆行列を計算する第３計算部と；前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する第４の計算部と；前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する第５の計算部と；前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する第６の計算部と；を備える。

前記（１）の発明によれば、大規模な連立方程式の行列演算等の煩雑で長時間を要する計算を行うことなく、簡易な計算により、極めて迅速で確実に、弾塑性材料の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。
前記（２）乃至（３）の発明によれば、前記（１）の発明よりも更に確実に弾塑性材料の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。
前記（４）乃至（５）の発明によれば、前記（２）乃至（３）の発明よりも更に確実に弾塑性材料の成形時におけるスプリングバックの発生部位を特定し、正確な成形に供することが可能となる。
前記（６）乃至（８）の発明によれば、（１）乃至（５）に記載のシミュレーション方法の結果に基づいて、複雑形状の加工部品であってもプレス成形時におけるスプリングバックの発生部位を効率良くピンポイントに特定し、その部位に起因するスプリングバックを抑えることで、寸法精度に優れた弾塑性材料の成形方法を提供することが可能となる。

Claims (16)

1. 弾塑性材料の成形シミュレーション方法であって、
   有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する工程と；
   計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する工程と；
   を含むことを特徴とする成形シミュレーション方法。
2. 計算された前記全等価節点力ベクトルに基づいて、前記弾塑性材料の全領域又は特定の領域のうちで、前記等価節点力ベクトルの大きい部位をスプリングバックの発生原因部位として特定する工程を更に含むことを特徴とする請求項１に記載の成形シミュレーション方法。
3. 全体剛性行列の逆行列を計算する工程と；
   前記全等価節点力ベクトルを外力ベクトルとして、前記全体剛性行列の逆行列と前記全等価節点力ベクトルとを乗算し、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの寄与度を計算する工程と；
   を更に含むことを特徴とする請求項１に記載の成形シミュレーション方法。
4. 前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの前記寄与度を表示する工程を更に含むことを特徴とする請求項３に記載の成形シミュレーション方法。
5. 全体剛性行列の逆行列を計算する工程と；
   前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する工程と；
   前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する工程と；
   前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する工程と；
   を更に含むことを特徴とする請求項１に記載の成形シミュレーション方法。
6. 前記第１の変位量と前記第２の変位量との間の変化量を表示する工程を更に含むことを特徴とする請求項５に記載の成形シミュレーション方法。
7. 請求項１に記載のシミュレーション方法の結果に基づいて、前記弾塑性材料を成形する弾塑性材料の成形方法。
8. 前記シミュレーション方法に基づいてスプリングバック発生原因部位を特定する工程と；
   前記スプリングバック発生原因部位に予め、エンボスを成形する工程と；
   前記エンボスに対し、圧縮応力を与えるように塑性変形させる工程と；
   を含む、請求項７に記載の、弾塑性材料の成形方法。
9. 前記エンボス成形工程で成形する前記エンボスが、２個以上の略同一寸法の円形エンボスであることを特徴とする、請求項８に記載の、弾塑性材料の成形方法。
10. 前記エンボス成形工程が、金型に設けられた凹部に勘合する高さ寸法の異なる複数の着脱式の工具を用いてエンボス形状の高さ寸法を調整する工程を有することを特徴とする、請求項８に記載の、弾塑性材料の成形方法。
11. コンピュータにより、弾塑性材料の成形シミュレーションを行うに際して、
    有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する工程と；
    計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する工程と；
    をコンピュータに実行させるためのプログラム。
12. 計算された前記全等価節点力ベクトルに基づいて、前記弾塑性材料の全領域又は特定の領域のうちで、前記等価節点力ベクトルの大きい部位をスプリングバックの発生原因部位として特定する工程を更にコンピュータに実行させるための請求項１１に記載のプログラム。
13. 全体剛性行列の逆行列を計算する工程と；
    前記全等価節点力ベクトルを外力ベクトルとして、前記全体剛性行列の逆行列と前記全等価節点力ベクトルとを乗算し、前記弾塑性材料の特定位置の変位量に対する前記外力ベクトルの成分ごとの寄与度を計算する工程と；
    を更にコンピュータに実行させるための請求項１１に記載のプログラム。
14. 全体剛性行列の逆行列を計算する工程と；
    前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する工程と；
    前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する工程と；
    前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する工程と；
    を更にコンピュータに実行させるための請求項１１に記載のプログラム。
15. 請求項１１に記載のプログラムを記録したコンピュータで読み取り可能な記録媒体。
16. 弾塑性材料の成形シミュレーションを行う成形シミュレーション装置であって、
    有限要素法を用いて、前記弾塑性材料の目標形状における１つ又は複数の有限要素ごとに、応力テンソルから要素等価節点力ベクトルを計算する第１の計算部と；
    計算された前記１つ又は複数の有限要素ごとの前記要素等価節点力ベクトルを、前記弾塑性材料の全領域又は特定の領域に亘って積分して、その領域の全等価節点力ベクトルを計算する第２の計算部と；
    全体剛性行列の逆行列を計算する第３計算部と；
    前記全等価節点力ベクトルを第１の外力ベクトルとして、前記全体剛性行列の逆行列と前記第１の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第１の変位量を計算する第４の計算部と；
    前記１つ又は複数の有限要素ごとに、前記第１の外力ベクトルからその要素の要素等価節点力ベクトルを除去した結果を第２の外力ベクトルとして、前記全体剛性行列の逆行列と前記第２の外力ベクトルとを乗算し、前記弾塑性材料の特定位置の第２の変位量を計算する第５の計算部と；
    前記１つ又は複数の有限要素ごとに、前記第１の変位量と前記第２の変位量との変化量を計算する第６の計算部と；
    を含むことを特徴とする成形シミュレーション装置。