【０００１】
【発明の属する技術分野】
本発明は、板状ブランクを用いたプレス成形の結果を有限要素法を用いてシミュレーションするプレス成形シミュレーションシステム、並びにプレス成形シミュレーション用プログラム、及びそのプログラムを記録した記録媒体に関し、コンピュータによる有限要素解析技術の分野に属する。
【０００２】
【従来の技術】
板材をブランクとし、これをプレス成形により所定形状の製品に成形する場合、まず、設計された製品形状を試作、検証した後、これを量産するための金型を設計、製作し、次いでこの金型を用いて製品を試作し、ＯＫであれば量産に移行する、というプロセスを経ることになるが、このようなプロセスにおいては、金型を用いた製品の試作の段階で、しわ、われ等の欠陥の発生が判明することがあり、この場合、製品の設計形状の見直しから金型の修正、試作まで、再度上記プロセスを実行しなければならないことになり、多大なコストと時間とを費やすことになる。
【０００３】
このような実情に対し、近年、製品設計の段階、或いは金型設計の段階で、成形結果をシミュレーションし、金型製造前のできるだけ早い段階で欠陥の発生を予測し、事前にその対策をとることにより上記のような多大な無駄を回避することが行われており、そのシミュレーションシステムとして各種のものが提供されている。
【０００４】
これらのシステムは、いずれも有限要素法を用いるものであって、非特許文献１によれば、その具体的解法として、動的陽解法、静的陽解法及び静的陰解法が挙げられている。
【０００５】
動的陽解法は、ブランクの変形に加速度を考慮するものであり、これに対して静的陽解法及び静的陰解法は、分割した各要素の変形を静的に捉え、各時間ステップごとに剛性マトリックスを組み立ててつりあい式を解いてゆくものであり、後者のものがプレス成形のシミュレーションには適しているとされている。
【０００６】
また、後者のうちの静的陽解法は、各時間ステップにおけるつりあい式を線形化して解くものであって、各ステップでの計算量は少なくてすむが、各ステップの増分を線形が維持されるように制限しなければならないため、精度のよいシミュレーションを行おうとするとステップ数が著しく多くなり、結果的に計算量が多くなるという欠点がある。
【０００７】
一方、静的陰解法は、変形を現実の変形現象に準じて非線形として捉え、各ステップごとに収束計算を行いながらステップを進めるものであって、高精度の解を得るために時間ステップの増分を小さく設定した微小増分法、計算効率の改善のために時間ステップの増分を大きく設定した大増分法があり、さらに、この静的陰解法に属するものとしてワンステップ法が存在する。
【０００８】
このワンステップ法は、製品の設計形状を平板状に展開したときのブランクの形状をワンステップで求めると共に、その際に各要素に生じる歪や応力から製品各部の板厚の分布やしわ、われ等の欠陥の発生を推測するものであり、計算量が少ないことや、金型に関する情報が不要であるなどの長所がある反面、成形中におけるブランクの変形過程を考慮しないので精度の点で課題がある。
【０００９】
【非特許文献１】
日本塑性加工学会編「わかりやすいプレス加工」日刊工業新聞社出版、
２０００年５月２３日発行、Ｐ１５９〜１７３
【００１０】
【発明が解決しようとする課題】
上記のように、プレス成形のシミュレーションシステムのための各解法はそれぞれ一長一短があり、例えば、製品の設計段階では、計算結果が速やかに得られるワンステップ法等を採用するシステムが用いられ、金型の設計、製作段階では、高精度な解が得られる静的陰解法の微小増分法を採用するシステムが用いられるが、前者の段階では、さらに精度を向上させることが要求され、後者の段階ではさらに計算速度の向上が要求されているところである。
【００１１】
そこで本発明は、有限要素法を用いたプレス成形のシミュレーションに関する以上のような実情に鑑み、この種のシミュレーションにおける精度と計算速度の両者を満足できる範囲で両立させ、製品の設計段階から金型の製作段階に至る広い場面で使用可能なシステムを提供することを課題とする。
【００１２】
【課題を解決するための手段】
上記課題を解決するため、本願各請求項の発明は次のように構成したことを特徴とする。
【００１３】
まず、請求項１に記載の発明は、板材をブランクとするプレス成形の結果を有限要素法を用いてシミュレーションするシステムに関するものであって、製品の設計データを入力する製品データ入力手段と、該入力手段で入力された製品データに基づき、製品形状を有限要素分割した解析モデルを作成する解析モデル作成手段と、該作成手段で作成された解析モデルを所定の平面に展開したときの各節点の該平面に直交する方向の直線変位量、及び各節点の法線ベクトルの向きを上記平面に直交する方向に一致させるのに要する回転変位量を計算する変位量計算手段と、該計算手段で計算された直線変位量及び回転変位量をそれぞれ複数個に分割する変位量分割手段と、上記解析モデルを平面形状へ展開したときの各要素の変形を、上記変位量分割手段で分割された複数の直線分割変位量及び回転分割変位量を各ステップの境界条件として複数ステップに分けて解析する解析手段と、該解析手段によって得られた平面形状での各要素の状態を製品形状に反映させるマッピング手段とを有することを特徴とする。
【００１４】
また、請求項２に記載の発明は、上記請求項１に記載のシステムにおいて、解析手段によって行われる解析計算のステップ数を設定するステップ数設定手段が設けられていることを特徴とする。
【００１５】
また、請求項３に記載の発明は、上記請求項１又は請求項２に記載のシステムにおいて、変位量分割手段は、変位量計算手段によって計算された各節点の直線変位量及び回転変位量を等分割することにより、複数の直線分割変位量及び回転分割変位量を生成することを特徴とする。
【００１６】
さらに、請求項４に記載の発明は、同じく請求項１又は請求項２に記載のシステムにおいて、各節点の直線変位及び回転変位の時間経過に対する特性を設定する変位特性設定手段が設けられ、変位量分割手段は、変位量計算手段によって計算された各節点の直線変位量及び回転変位量を上記変位特性設定手段で設定された変位特性に基づいて分割することにより、複数の直線分割変位量及び回転分割変位量を生成することを特徴とする。
【００１７】
そして、請求項５に記載の発明は、上記請求項１から請求項４のいずれかに記載のシステムにおいて、解析手段は、各ステップの解析を静的陰解法によって行うことを特徴とする。
【００１８】
以上の請求項１から請求項５に記載のシステムによれば、変位量計算手段によって計算された解析モデルを平面形状に展開したときの各節点の直線変位量及び回転変位量が変位量分割手段によってそれぞれ複数の分割変位量に分割されると共に、解析手段により、これらの分割変位量をそれぞれの境界条件とする複数のステップで、解析モデルを平面形状に展開したときの変形の解析計算が行われることになる。
【００１９】
これは、解析モデルを所定の境界条件のもとでワンステップで展開形状へ変形させる前述のワンステップ法を、例えば１０回程度の複数回に分けて実行することに相当し、したがって、ワンステップ法の長所である計算効率のよさを生かしながら、その解析精度が向上することになる。
【００２０】
その場合に、請求項２に記載のシステムによれば、解析手段にって行う解析計算のステップ数をステップ数設定手段によって任意に設定することができるから、製品形状に応じて或いは精度と計算速度のいずれを重視するか等に応じて、最適なステップ数を設定することが可能となり、精度と計算速度の一層良好な両立が可能となる。
【００２１】
また、請求項３に記載のシステムによれば、変位量分割手段によって直線分割変位量及び回転分割変位量を生成する場合に、変位量計算手段で計算された各節点の直線変位量及び回転変位量を等分割するだけであるから、複雑な処理を要することなく、各分割変位量を容易に生成することができ、それだけシステムの構成や処理内容が簡素化されることになる。
【００２２】
これに対して、請求項４に記載のシステムによれば、予め各節点の直線変位及び回転変位の時間経過に対する特性を設定しておき、変位量分割手段によって直線分割変位量及び回転分割変位量を生成する場合に、各節点の直線変位量及び回転変位量を上記変位特性設定手段で設定された変位特性に基づいて分割するので、各ステップの解析計算が各節点ごとに適切に設定された境界条件のもとで行われることになり、したがって、最終的な解析結果の精度がさらに向上することになる。
【００２３】
そして、請求項５に記載のシステムによれば、各ステップの解析計算が静的陰解法によって行われ、良好な解析精度が得られる。
【００２４】
一方、請求項６に記載の発明は、板材をブランクとするプレス成形の結果を有限要素法を用いてシミュレーションするためのプログラムに関するものであって、与えられた製品形状を有限要素分割して解析モデルを作成する解析モデル作成機能、境界条件を取得して解析モデルの各要素の変形を解析計算する解析機能、及び平面形状での各要素の状態を製品形状に反映させるマッピング機能を備えたコンピュータで用いられ、該コンピュータを、上記解析モデル作成機能によって作成された解析モデルを所定の平面に展開したときの各節点の該平面に直交する方向の直線変位量及び各節点の法線ベクトルの向きを上記平面に直交する方向に一致させるのに要する回転変位量を計算する変位量計算手段、並びに上記解析機能による解析計算を複数ステップに分けて行う場合の各ステップの境界条件とするために、上記変位量計算手段で計算された直線変位量及び回転変位量を複数個に分割する変位量分割手段として機能させることを特徴とする。
【００２５】
また、請求項７に記載の発明は、上記請求項６に記載のプログラムにおいて、コンピュータを、解析計算のステップ数を設定するステップ数設定手段として機能させると共に、変位量分割手段として機能させるときに、直線変位量及び回転変位量を上記ステップ数設定手段としての機能によって設定されたステップ数と同数個に分割することを特徴とする。
【００２６】
また、請求項８に記載の発明は、上記請求項６または請求項７に記載のプログラムにおいて、コンピュータを、変位量分割手段として機能させるときに、変位量計算手段によって計算された各節点の直線変位量及び回転変位量を等分割することを特徴とする。
【００２７】
また、請求項９に記載の発明は、同じく請求項６または請求項７に記載のプログラムにおいて、コンピュータを、各節点の直線変位及び回転変位の時間経過に対する特性を設定する変位特性設定手段として機能させると共に、変位量分割手段として機能させるときに、変位量計算手段によって計算された各節点の直線変位量及び回転変位量を上記変位特性設定手段で設定された変位特性に基づいて分割することを特徴とする。
【００２８】
以上の請求項６から請求項９に記載の発明は、上記請求項１から請求項４に記載のシステムに関する発明に対応するもので、当該プログラムを、解析モデル作成機能、解析機能、及びマッピング機能等を備えたコンピュータに適用することにより、請求項１から請求項４に記載のシステムと同様の作用が得られる。
【００２９】
さらに、請求項１０に記載の発明は、請求項６から請求項９のいずれかに記載のプログラムを記録したコンピュータ読み取り可能な記録媒体に関し、この発明によれば、当該記録媒体を、解析モデル作成機能、解析計算機能、及びマッピング機能等を備えたコンピュータに適用することにより、請求項１から請求項４に記載のシステムと同様の作用が得られる。
【００３０】
【発明の実施の形態】
次に、本発明の実施の形態に係るプレス成形のシミュレーションシステムについて説明する。
【００３１】
図１は、このシステムの中心となるコンピュータの構成を示すもので、このコンピュータ１０は、中央処理装置１１と、各種条件の設定やシステムの制御等に用いられる入力装置１２と、記録媒体２０からのデータ読込み装置１３と、プログラム記録部１４ａ及びデータファイル記録部１４ｂが設けられた記録装置１４と、入力画面や計算結果等を表示する表示装置１５と、計算結果等を印刷する印刷装置１６とを有する。
【００３２】
上記記録装置１４のプログラム記録部１４ａには、図２に示すように、メインプログラムと複数のサブプログラムとが記録され、サブプログラムとしては、解析モデル作成サブプログラム、基礎データ作成サブプログラム、解析計算サブプログラム、マッピング処理サブプログラム及び描画サブプログラムが記録されている。
【００３３】
また、上記記録装置１４のデータファイル記録部１４ｂには、図３に示すように、解析モデルデータファイルと、解析用基礎データファイルと、材料属性データファイルと、解析結果データファイルとが記録されるようになっている。
【００３４】
これらのデータファイルのうち、解析モデルデータファイルは、解析対象製品を有限要素分割してなる解析モデルについての各種のデータを記録したテーブルで構成され、図４に示す各節点の座標データを記録した節点座標テーブル、図５に示す各要素を構成する節点を記録した要素構成テーブル、及び図６に示す各節点を構成する要素を記録した節点構成テーブルとを含んでいる。
【００３５】
また、解析用基礎データファイルは、解析計算に先立ち、予め計算され或いは設定された各種のデータを記録したテーブルで構成され、図７に示す要素ベクトルテーブル、図８に示す節点ベクトルテーブル、図９に示す直線変位量テーブル、図１０に示す回転変位量テーブル、図１１に示す拘束条件テーブル、図１２に示す変位特性曲線テーブル等を含んでいる。
【００３６】
また、材料属性データファイルは、図１３に示す解析計算に必要な解析対象製品の各種の材料属性を記録したテーブルによって構成されている。
【００３７】
さらに、解析結果データファイルは、図１４に示すように、複数の解析ステップごとに解析モデルの全要素についての応力、歪及び厚さ等の解析結果データを記録するテーブルによって構成されている。
【００３８】
次に、当該システムの作用を、上記記録装置１４のプログラム記録部１４ａに記録されているメインプログラムの動作を示す図１５のフローチャートに従って説明する。
【００３９】
まず、プロセスＰ１として、図１に示すコンピュータ１０のデータ読み込み装置１３により、記録媒体２０から解析対象の製品のＣＡＤデータを読み込む。このとき、コンピュータ１０の表示装置１５には、図１６に示すように、当該製品の全体形状を表示した画面１５ａが表示される。
【００４０】
次に、プロセスＰ２として、上記プログラム記録部１４ａに記録されている解析モデル作成サブプログラムにより、上記ＣＡＤデータに基づいて当該製品の形状が有限要素分割されてなる解析モデルが作成される。このとき、図１７に示すように、作成された解析モデルが表示装置１５に画面１５ｂとして表示されると共に、図１８に示すように、解析モデルを構成するメッシュの各節点に節点番号Ｎ１、Ｎ２…が付され、また各要素に要素番号Ｅ１、Ｅ２…が付される。
【００４１】
そして、上記解析モデル作成サブプログラムにより、解析モデルを所定のｘｙｚ座標系に置いたときの各節点Ｎ１、Ｎ２…のｘ座標、ｙ座標、ｚ座標が読み取られ、図４の節点座標テーブルに各節点ごとに書き込まれる。また、各要素Ｅ１、Ｅ２…を構成する節点、即ち当該要素の周囲の節点の番号が読み取られ、図５に示す要素構成テーブルに各要素ごとに書き込まれ、さらに、各節点Ｎ１、Ｎ２…を構成する要素、即ち当該節点を取り囲む要素の番号が読み取られ、図６の節点構成テーブルに各節点ごとに書き込まれる。
【００４２】
次に、プロセスＰ３として、上記解析モデルについての解析用基礎データが上記プログラム記録部１４ａに記録されている基礎データ作成サブプログラムにより作成される。
【００４３】
ここで、このサブプログラムの動作を図１９に示すフローチャートに従って説明する。
【００４４】
まず、プロセスＰ２１で、前述の解析モデル作成処理により得られた節点の座標データ、要素構成データ及び節点構成データを図４、図５、図６のテーブルからそれぞれ読み出し、次に、プロセスＰ２２で、図１８に示す各要素Ｅ１、Ｅ２…についての法線ベクトル（以下、「要素ベクトル」という）ＶＥ１、ＶＥ２…を算出する。
【００４５】
この要素法線ベクトルＶＥ１、ＶＥ２…の算出方法を例えば要素Ｅ１について説明すると、まず図５の要素構成テーブルから要素Ｅ１が節点Ｎ１、Ｎ２、Ｎ５、Ｎ４で構成されることを読み出すと共に、これらの節点Ｎ１、Ｎ２、Ｎ５、Ｎ４のｘ座標、ｙ座標、ｚ座標を図４の節点座標テーブルから読み出し、これらの座標データに基づき、要素Ｅ１を平面としたときの該要素Ｅ１の傾きを求める。そして、その平面に直行する方向の単位ベクトルのｘ成分、ｙ成分、ｚ成分を求め、これらの成分を有するベクトルを要素Ｅ１の要素ベクトルＶＥ１とする。
【００４６】
このようにして各要素Ｅ１、Ｅ２…についての要素ベクトルＶＥ１、ＶＥ２…のｘ成分、ｙ成分、ｚ成分を順次求め、これを図７の要素ベクトルテーブルに書き込む。
【００４７】
次に、プロセスＰ２３で、各節点Ｎ１、Ｎ２…についての法線ベクトル（以下、「節点ベクトル」という）ＶＮ１、ＶＮ２…を算出する。
【００４８】
この節点法線ベクトルＶＮ１、ＶＮ２…の算出方法を例えば図１８に示す節点Ｎ５について説明すると、図６の節点構成テーブルから節点Ｎ５を構成する要素、即ち節点Ｎ５を取り囲む要素を読み出す。この場合、図１８の例では、要素Ｅ１、Ｅ２、Ｅ４、Ｅ３が読み出される。
【００４９】
そして、図７の要素ベクトルテーブルから、これらの要素Ｅ１、Ｅ２、Ｅ４、Ｅ３についての要素ベクトルＶＥ１、ＶＥ２、ＶＥ４、ＶＥ３のｘ成分、ｙ成分、ｚ成分をそれぞれ読み出すと共に、これら４つのベクトルの各成分の平均値をｘ成分、ｙ成分、ｚ成分とする単位ベクトルを求め、これを節点Ｎ５についての節点ベクトルＶＮ５とする。
【００５０】
このようにして、各節点Ｎ１、Ｎ２…についての節点ベクトルＶＮ１、ＶＮ２…を順次求め、そのｘ成分、ｙ成分、ｚ成分を図８の節点ベクトルテーブルに書き込む。
【００５１】
次に、プロセスＰ２４で、全節点Ｎ１、Ｎ２…についての節点ベクトルＶＮ１、ＶＮ２…から、これらのベクトルの方向の平均的な方向として、解析モデルの展開方向、即ち解析モデルを押しつぶして平面形状とする際の押圧方向を決定し、図２０に示すように、その方向を例えばＺ方向とするＸＹＺ座標系を設定する。そして、このＸＹＺ座標系上で展開方向Ｚに直交する平面を設定し、これを解析モデルの仮想展開面Ｐをする。その場合に、この展開面ＰのＺ方向の位置は任意に設定することがでる。
【００５２】
次に、プロセスＰ２５で、解析モデルを上記仮想展開面Ｐ上へ展開したときの各節点Ｎ１、Ｎ２…の展開方向Ｚの直線変位量Ｄ１、Ｄ２…を算出し、これを図９に示す直線変位量テーブルに展開方向と共に書き込む。
【００５３】
その場合に、この直線変位量Ｄ１、Ｄ２…の計算は、図４の節点座標テーブルに記録されている各節点Ｎ１、Ｎ２…のｘｙｚ座標系での座標を上記展開方向をＺ方向とするＸＹＺ座標系での座標に変換した上で、このＸＹＺ座標系で、各節点Ｎ１、Ｎ２…から仮想展開面ＰまでのＺ方向の距離を求めることにより行われ、図２０に節点Ｎ１についての直線変位量Ｄ１を図示している。
【００５４】
また、プロセスＰ２６で、解析モデルを上記仮想展開面Ｐ上へ展開したときの各節点ベクトルＶＮ１、ＶＮ２…の軸Ｘ、Ｙ、Ｚ周りの各回転変位量（ＡＸ１、ＡＹ１、ＡＺ１）、（ＡＸ２、ＡＹ２、ＡＺ２）…を算出し、これを図１０の回転変位量テーブルに書き込む。
【００５５】
この回転変位量は、節点ベクトルの向きを解析モデルの展開方向に一致させるのに要する各軸周りの回転角であって、図２０により、節点Ｎ１の節点ベクトルＶＮ１を例にとって説明すると、このベクトルＶＮ１をＹＺ平面に平行な面上に移動させるのに要するＹ軸周りに回転角をＡＹ１とし、また、このベクトルをさらにＺＸ平面に平行な面上に移動させるのに要するＸ軸周りの回転角をＡＸ１とすれば、これらの回転角により、節点ベクトルＶＮ１の向きはＺ方向に一致することになる。
【００５６】
つまり、これらの回転角ＡＹ１、ＡＸ１がベクトルＶＮ１を展開方向に一致させるためのＹ軸及びＸ軸周りの回転変位量となる。なお、この場合、展開方向とＺ方向とが一致しているから、Ｚ軸周りの回転変位量は０となる。
【００５７】
次に、プロセスＰ２７で、各節点Ｎ１、Ｎ２…について、直線変位と回転変位の変位特性曲線を定義し、その曲線を示すデータを図１２の変位特性曲線テーブルに記録する。
【００５８】
この変位特性曲線は、解析モデルを仮想展開面Ｐに展開する際の直線変位と回転変位の時間経過に対する変化の特性を示すものであり、例えば次のように設定される。
【００５９】
図２１に示すように、今、例えば解析モデルの節点Ｎａが節点Ｎｂに対して符号アで示す位置にあり、これが平面形状に展開されたときに符号イで示す位置に変位するものとする。このときの節点Ｎａの最終的な直線変位量をＤ、節点ベクトルＶＮａの紙面に直交する軸周りの回転変位量をＡとする。また、符号ウで示す変形過程の中間の段階までの節点Ｎａの直線変位量をｄ、節点ベクトルＶＮａの回転変位量をａとし、節点Ｎａ、Ｎｂ間の距離をＲとすれば、変形過程での直線変位量ｄと回転変位量ａの関係は、
ｄ＝Ｄ−Ｒ・ｓｉｎ（Ａ−ａ）
Ｒ＝Ｄ／ｓｉｎＡ
であるから、
ｄ＝Ｄ・［１−ｓｉｎ（Ａ−ａ）／ｓｉｎＡ］
となる。
【００６０】
ここで、例えば最終直線変位量Ｄを３、最終回転変位量Ａを０．７８５（４５°）とすると共に、中間段階の回転変位量ａとして、最終回転変位量Ａを１０等分した値、０．０７９、０．１５７…を設定し、これらの値についての各段階の直線変位量を計算すれば、図２２の表に示す値が得られる。
【００６１】
また、例えば時間経過に対して回転変位が直線的に推移するものとして、上記の計算結果を時間経過に対する特性として示せば、図２３に示すような直線変位特性曲線ＬＤ、及び回転変位特性曲線ＬＡが得られる。
【００６２】
ここで、実際には時間経過に対する直線変位とＸ軸周りの回転変位の特性とＹ軸周りの回転変位の特性とを定義しておく必要があり、その場合に、直線変位が時間経過に対して直線的に推移するものとし、そのときのＸ軸周り及びＹ軸周りの回転変位がどのような特性で変化するかを計算して定義しておくようにしてもよい。
【００６３】
また、最終変位量が相違しても時間経過に対して同様の特性で変形する節点が多数存在するので、変位特性を定義する際には、最終変位量を１とし、各時間ごとの中間段階の変位量については全変位量に対する比率で定義して、これを特定節点の変位特性として使用するときに、特性曲線上の値（比率）にその節点の最終変位量を掛けることにより各段階の中間変位量を求めるようにしてもよい。このようにすれば、変位特性が共通の複数の節点について変位特性曲線が共用され、その設定作業や記録容量が節減されることになる。
【００６４】
具体的には、図１１の変位特性曲線テーブルに示すように、特性曲線番号をキーとして、解析時間を１０等分したときの各時間までの直線変位、Ｘ軸周りの回転変位及びＹ軸周りの回転変位量を、最終変位量を１としたときの比率で記録すると共に、図９、図１０の各節点ごとの直線変位量及び回転変位量をそれぞれ記録した直線変位量テーブル及び回転変位量テーブルに、それぞれの節点について適用される変位特性曲線の番号を記録する。
【００６５】
次に、プロセスＰ２８、Ｐ２９で、解析計算時の拘束条件として、解析モデルの展開時にＸ方向及びＹ方向への変位をそれぞれ禁止した２つの節点を指定する。これは、解析モデルを平面形状に展開するときに、すべての節点がＸ方向及びＹ方向に移動可能であると、ＸＹ平面上でのモデル全体としての位置が決まらず、計算不能となることを防止するためである。
【００６６】
そして、例えば節点ＮｘをＸ方向への変位を禁止した節点、節点ＮｙをＹ方向への変位を禁止した節点に指定したものとすれば、図１１の拘束条件テーブルに、それらの節点番号Ｎｘ、Ｎｙと、その節点のＸ方向及びＹ方向の変位がそれぞれ０であることを記録する。
【００６７】
以上により、基礎データ作成サブプログラムによる基礎データの作成処理が終了し、プロセスＰ３０で、得られたデータが解析用基礎データファイルに格納される。そして、コンピュータ１０の動作は、図１５にフローチャートを示すメインプログラムによる処理に戻り、プロセスＰ４以下を実行する。
【００６８】
メインプログラムのプロセスＰ４では、当該コンピュータ１０の入力装置１２からの入力に基づいて解析計算のステップ数を設定する。また、プロセスＰ５では、今回の解析対象の製品の材料に関するデータ、即ち当該材料のヤング率、ポアソン比、降伏応力、接線係数等の解析に必要なデータを図１３のテーブルから読み取る。そして、プロセスＰ６〜Ｓ１０の解析処理を開始する。
【００６９】
この解析処理においては、まずプロセスＰ６で、第１ステップの解析計算に必要な境界条件を作成する。
【００７０】
つまり、図９、図１０のテーブルから、解析モデルを平面形状まで展開したときの各節点Ｎ１、Ｎ２…の最終的な直線変位量Ｄ１、Ｄ２…、各節点ベクトルＶＮ１、ＶＮ２…のＸ軸周り及びＹ軸周りの最終的な回転変位量ＡＸ１、ＡＸ２…、ＡＹ１、ＡＹ２…、並びにこれらの変位の変位特性曲線番号ＬＤ１…、ＬＡＸ１…、ＬＡＹ１…を読み出す。
【００７１】
そして、上記プロセスＰ５で設定された解析ステップ数が例えば１０であるとすれば、図１２のテーブルで上記番号ＬＤ１…、ＬＡＸ１…、ＬＡＹ１…をキーとして定義された各変位特性曲線の解析時間を１０等分したときの最初の時間での値（比率）を読み出し、その値と、上記直線変位量Ｄ１、Ｄ２…、Ｘ軸周り及びＹ軸周りの回転変位量ＡＸ１、ＡＸ２…、ＡＹ１、ＡＹ２…とを、各節点ごとに掛けることにより、各節点の第１ステップでの境界条件を算出する。
【００７２】
これにより、各節点Ｎ１、Ｎ２の直線変位の第１ステップでの境界条件、各節点ベクトルＶＮ１、ＶＮ２…のＸ軸周り及びＹ軸周りの回転変位の第１ステップでの境界条件が得られることになる。
【００７３】
ここで、上記の例は、各変位の特性曲線が１０個のデータで定義され、かつ設定された解析ステップ数も１０であるから、特性曲線の定義で用いられたデータをそのまま境界条件の算出に用いることができるが、これらが一致しない場合は、特性曲線から補間法を用いて各ステップの境界条件を算出することになる。
【００７４】
このようにしてすべての節点についての直線変位、Ｘ軸周りの回転変位及びＹ軸周りの回転変位の境界条件が得られると、次にメインプログラムのプロセスＰ７で、図４に示す各節点Ｎ１、Ｎ２…の解析モデル上でのｘ座標、ｙ座標、ｚ座標をＸＹＺ座標系に変換した値に基づいて初期条件を設定し、これに上記のようにして得られた境界条件と、図１１のテーブルに示す２つの節点Ｎｘ、Ｎｙについての拘束条件とを用い、解析計算サブプログラムにより、第１ステップの解析計算を実行する。
【００７５】
この計算は静的陰解法によって行われ、ヤング率やポアソン比等の材料属性データを用いて剛性マトリックスを設定し、各要素ごとの釣り合い式から、各要素Ｅ１、Ｅ２…の応力、歪、厚さ等の変化を計算する。そして、プロセスＰ８で、その計算結果を図１４のテーブルに第１ステップの全要素Ｅ１、Ｅ２…についての解析結果として記録する。
【００７６】
次に、プロセスＰ９で、全ステップの解析計算が終了したか否かを判定し、現時点では未だ終了していないから、次にプロセスＰ１０を実行し、プロセスＰ８で解析結果テーブルに記録した第１ステップの解析結果を読み出し、これに基づいて第２ステップの解析計算の初期条件に設定すると共に、プロセスＰ６を再び実行し、各節点についての直線変位、Ｘ軸周り及びＹ軸周りの回転変位の境界条件を作成する。
【００７７】
この場合は、図９、図１０のテーブルから読み取った各節点Ｎ１、Ｎ２…の最終的な直線変位量Ｄ１、Ｄ２…、各節点ベクトルＶＮ１、ＶＮ２…のＸ軸周り及びＹ軸周りの最終的な回転変位量ＡＸ１、ＡＸ２…、ＡＹ１、ＡＹ２…と、図１２のテーブルで定義された各変位特性曲線の解析時間を１０等分したときの２番目の時間での値（比率）とを掛けることにより、各節点ないし節点ベクトルの第２ステップでの境界条件を求める。
【００７８】
そして、これらの境界条件と、プロセスＰ１０で設定した前回の解析結果に基づいて設定した初期条件と、図１１のテーブルに記録されている拘束条件とを用い、プロセスＰ７で、解析計算サブプログラムにより再び解析計算を実行し、プロセスＰ８で、その計算結果を図１４の解析結果テーブルに第２ステップの全要素Ｅ１、Ｅ２…についての解析結果として記録する。
【００７９】
以下同様にして、プロセスＰ９で全ステップの解析計算が終了したものと判定するまでプロセスＰ６〜Ｐ１０を繰り返し、前回の解析結果に基づいて設定した初期条件と、変位特性曲線の解析時間をステップ数で分割した各時間での値を用いて算出した境界条件と、テーブルに記録されている拘束条件とを用いて解析計算を繰り返し実行すると共に、得られた計算結果を各ステップごとに図１４の解析結果テーブルに記録する。
【００８０】
このようにして、設定されたステップ数の解析計算が終了すれば、解析モデルは平面形状に展開されたことになり、次に、プロセスＰ９からプロセスＰ１１を実行し、マッピング処理を行う。
【００８１】
この処理は、マッピング処理サブプログラムによって実行され、図１４の解析結果テーブルに記録されている最終ステップの計算結果、即ち平面形状に展開した状態での各要素Ｅ１、Ｅ２…の応力、歪、厚さ等のデータを、展開前の解析モデルの各要素Ｅ１、Ｅ２…に当てはめることにより、一様な板状ブランクを当該製品形状にプレス成形した際の各部の応力、歪、厚さ等の分布を求めるものである。
【００８２】
そして、プロセスＰ１２で、描画サブプログラムにより、例えば図２４の画面１５ｃに示すように、マッピング処理によって得られた製品形状上での厚さの分布等を表示する。
【００８３】
以上のようにして、板材からＣＡＤデータによって与えられた製品形状へのプレス成形の状況がシミュレーションされることになり、しわ、われ等のプレス成形時の欠陥を事前に推測することが可能となる。
【００８４】
なお、以上の実施の形態では、各節点の直線変位と、Ｘ軸周り及びＹ軸周りの回転変位の時間経過に対する特性を定義し、これを分割して各ステップの境界条件を求めるようにしたが、システム構成の簡素化や計算の迅速性が要求される場合には、直線変位及び回転変位のいずれもが時間経過に対して直線的に推移するものとして処理するようにしてもよく、この場合、各ステップでの境界条件は、最終的な変位量を設定されたステップ数で等分割したときの値となる。
【００８５】
また、以上の説明のように、この実施の形態に係るシステムは複数のプログラムによって動作するように構成されているが、そのうち、有限要素分割プログラム、解析計算プログラム、マッピング処理プログラム、描画プログラム等は既存のものを利用することができる。したがって、これらの既存プログラムがインストールされている汎用のシステムに、基礎データ作成プログラムと、これらのプログラムを関連付けるためのプログラムを読み込ませるだけで、当該システムを本発明を実施するシステムとして動作させることが可能となり、基礎データ作成プログラムにより、予め製品ごとに基礎データを作成しておくことも可能である。
【００８６】
その場合に、新たに読み込ませる基礎データ作成プログラム等は、本発明に係るシミュレーション用プログラムの実施の形態を構成することになり、これを記録したコンピュータ読み取り可能な記録媒体は、本発明に係る記録媒体の実施の形態を構成することになる。
【００８７】
【発明の効果】
以上のように、本発明に係るシステム、プログラム及び記録媒体によれば、有限要素法を用いたプレス成形のシミュレーションシステムとして、従来のワンステップ法の特長である計算速度の速さに加えて、シミュレーションの精度が向上する。したがって、製品の設計段階から金型の製作段階に至る広い場面で使用可能なシステムが実現されることになる。
【００８８】
また、製品形状を有限要素分割する解析モデル作成プログラムや解析計算プログラム等は既存のものを利用することができるので、既存の汎用システム等を利用して上記効果を実現することができる利点がある。
【図面の簡単な説明】
【図１】本発明の実施の形態に係るシステムの全体構成を示すブロック図である。
【図２】同システムの記録装置に記録されるプログラムの説明図である。
【図３】同じくデータファイルの説明図である。
【図４】有限要素データファイルに含まれる節点座標テーブルの説明図である。
【図５】同じく要素構成テーブルの説明図である。
【図６】同じく節点構成テーブルの説明図である。
【図７】解析用基礎データファイルに含まれる要素ベクトルテーブルの説明図である。
【図８】同じく節点ベクトルテーブルの説明図である。
【図９】同じく直線変位量テーブルの説明図である。
【図１０】同じく回転変位量テーブルの説明図である。
【図１１】同じく変位特性曲線テーブルの説明図である。
【図１２】同じく拘束条件テーブルの説明図である。
【図１３】材料属性データファイルを構成するテーブルの説明図である。
【図１４】解析結果データファイルを構成するテーブルの説明図である。
【図１５】システム全体を制御するメインプログラムの動作を示すフローチャートである。
【図１６】ＣＡＤ図面を表示する画面の説明図である。
【図１７】有限要素分割した解析モデルを表示する画面の説明図である。
【図１８】解析モデルの要素や節点等の説明図である。
【図１９】解析用基礎データ作成プログラムの動作を示すフローチャートである。
【図２０】解析モデル展開時の変位量の説明図である。
【図２１】変位特性曲線の定義方法の一例を示す説明図である。
【図２２】その方法で得られた変位特性曲線のデータを示すテーブルである。
【図２３】その変位特性曲線を図示したグラフである。
【図２４】シミュレーション結果を示す画面の説明図である。
【符号の説明】
１０ コンピュータ
１１ 中央処理装置
１２ 入力装置
１３ データ読み込み装置
１４ 記憶装置
１５ 表示装置[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a press forming simulation system that simulates the result of press forming using a plate-shaped blank using a finite element method, a press forming simulation program, and a recording medium on which the program is recorded, and a finite element analysis by a computer. Belongs to the field of technology.
[0002]
[Prior art]
When a plate material is used as a blank and formed into a product of a predetermined shape by press molding, first, a prototype of the designed product shape is verified and verified, and then a mold for mass-producing the product is designed and manufactured. The process of trial production of a product using a mold and shifting to mass production if it is OK goes through. In such a process, wrinkles, we, etc., occur at the stage of trial production of a product using a mold. In some cases, the occurrence of defects may be found. In this case, the above process must be performed again from the review of the product's design shape to the correction of the mold and the trial production, which consumes a great deal of cost and time. Will be.
[0003]
In order to cope with such a situation, in recent years, a molding result is simulated at a product design stage or a mold design stage, a defect occurrence is predicted as early as possible before mold production, and a countermeasure is taken in advance. As a result, a large amount of waste as described above is avoided, and various simulation systems have been provided.
[0004]
Each of these systems uses the finite element method, and according to Non-Patent Document 1, as a specific solution, a dynamic explicit method, a static explicit method, and a static implicit method are mentioned.
[0005]
The dynamic explicit method considers the acceleration of the deformation of the blank, whereas the static explicit method and the static implicit method statically capture the deformation of each divided element and provide a rigid matrix at each time step. Are assembled to solve the balance equation, and the latter is said to be suitable for press molding simulation.
[0006]
The static explicit method of the latter solves the equilibrium equation at each time step by linearizing it, so that the amount of calculation at each step is small, but the linearity is maintained in the increment of each step. Therefore, if an attempt is made to perform an accurate simulation, the number of steps is significantly increased, resulting in a disadvantage that the amount of calculation is increased.
[0007]
On the other hand, in the static implicit method, the deformation is treated as nonlinear according to the actual deformation phenomenon, and the steps are advanced while performing the convergence calculation for each step.In order to obtain a highly accurate solution, the time step is increased. There is a small increment method in which is set small, a large increment method in which the increment of the time step is set large to improve the computational efficiency, and a one-step method belongs to the static implicit method.
[0008]
In this one-step method, the shape of the blank when the design shape of the product is developed into a flat plate is determined in one step, and the distribution and wrinkles of the plate thickness of each part of the product are determined from the strain and stress generated in each element at that time. It is an estimate of the occurrence of defects such as defects, and has advantages such as a small amount of calculation and no need for information on the mold.However, since it does not consider the blank deformation process during molding, there is a problem in accuracy. There is.
[0009]
[Non-patent document 1]
The Japan Society for Technology of Plasticity, “Easy Pressing”, published by Nikkan Kogyo Shimbun,
Published on May 23, 2000, P159-173
[0010]
[Problems to be solved by the invention]
As described above, each solution for the press molding simulation system has advantages and disadvantages. For example, in the product design stage, a system that employs a one-step method or the like in which calculation results are quickly obtained is used, and a mold is used. In the design and production stages, a system that employs the small incremental method of the static implicit method that can obtain a highly accurate solution is used, but in the former stage, it is required to further improve the accuracy, and in the latter stage, Further improvement in calculation speed is being demanded.
[0011]
Accordingly, the present invention has been made in view of the above-described circumstances regarding the simulation of press forming using the finite element method, and has achieved both the accuracy and the calculation speed in this type of simulation within a range that can be satisfied. It is an object of the present invention to provide a system that can be used in a wide range of situations from the production stage to the production stage.
[0012]
[Means for Solving the Problems]
In order to solve the above problems, the invention of each claim of the present application is characterized in that it is configured as follows.
[0013]
First, the invention according to claim 1 relates to a system for simulating the result of press forming using a plate material as a blank using a finite element method, and a product data input means for inputting product design data, Analysis model creating means for creating an analytic model obtained by dividing the product shape into finite elements based on the product data input by the input means; and an analysis model created by the creating means, the nodes of each of the nodes when the analysis model is developed on a predetermined plane. Displacement amount calculating means for calculating a linear displacement amount in a direction perpendicular to the plane and a rotational displacement amount required to make the direction of a normal vector of each node coincide with the direction perpendicular to the plane; Displacement amount dividing means for dividing the linear displacement amount and the rotational displacement amount into a plurality of pieces respectively, and the deformation of each element when the analysis model is developed into a planar shape, the displacement amount Analyzing means for analyzing a plurality of linearly divided displacement amounts and rotationally divided displacement amounts divided by the dividing means into a plurality of steps as boundary conditions for each step, and a state of each element in a planar shape obtained by the analyzing means And a mapping means for reflecting the symbol on the product shape.
[0014]
According to a second aspect of the present invention, in the system of the first aspect, a step number setting means for setting the number of steps of the analysis calculation performed by the analysis means is provided.
[0015]
According to a third aspect of the present invention, in the system according to the first or second aspect, the displacement amount dividing means calculates a linear displacement amount and a rotational displacement amount of each node calculated by the displacement amount calculating means. By dividing equally, a plurality of linearly divided displacement amounts and rotationally divided displacement amounts are generated.
[0016]
Further, the invention according to claim 4 is the system according to claim 1 or 2, further comprising a displacement characteristic setting means for setting a characteristic of each node with respect to a lapse of time of a linear displacement and a rotational displacement. The amount dividing unit divides the linear displacement amount and the rotational displacement amount of each node calculated by the displacement amount calculating unit based on the displacement characteristic set by the displacement characteristic setting unit, so that a plurality of linear divided displacement amounts and It is characterized in that a rotation split displacement amount is generated.
[0017]
According to a fifth aspect of the present invention, in the system according to any one of the first to fourth aspects, the analysis means performs an analysis of each step by a static implicit method.
[0018]
According to the above-described system, the linear displacement amount and the rotational displacement amount of each node when the analytical model calculated by the displacement amount calculating means is developed into a planar shape are calculated by the displacement amount dividing means. And the analysis means performs an analysis calculation of the deformation when the analysis model is developed into a planar shape in a plurality of steps using these divided displacement amounts as respective boundary conditions. Will be
[0019]
This is equivalent to executing the above-described one-step method of deforming the analysis model into a developed shape in one step under predetermined boundary conditions in a plurality of times, for example, about 10 times. The analysis accuracy will be improved while taking advantage of the computational efficiency, which is an advantage of the method.
[0020]
In that case, according to the system of claim 2, the number of steps of the analysis calculation performed by the analysis means can be arbitrarily set by the number-of-steps setting means. The optimum number of steps can be set according to which of the speeds is to be emphasized, and the accuracy and the calculation speed can be both better balanced.
[0021]
According to the system of claim 3, when the displacement dividing means generates the linear divided displacement and the rotational divided displacement, the linear displacement and the rotational displacement of each node calculated by the displacement computing means. Since the amount is simply divided equally, each divided displacement amount can be easily generated without requiring complicated processing, and the configuration and processing contents of the system are simplified accordingly.
[0022]
On the other hand, according to the system of the fourth aspect, the characteristics of the respective nodes with respect to the lapse of time of the linear displacement and the rotational displacement are set in advance, and the displacement dividing means sets the linear divided displacement and the rotational divided displacement. Is generated, the linear displacement amount and the rotational displacement amount of each node are divided based on the displacement characteristics set by the displacement characteristic setting means, so that the analytical calculation of each step is appropriately set for each node. It will be performed under boundary conditions, thus further improving the accuracy of the final analysis results.
[0023]
According to the system of the fifth aspect, the analytical calculation of each step is performed by the static implicit method, and good analytical accuracy can be obtained.
[0024]
On the other hand, the invention according to claim 6 relates to a program for simulating the result of press forming using a plate material as a blank by using a finite element method, and analyzes a given product shape by dividing it into finite elements. A computer equipped with an analysis model creation function for creating a model, an analysis function for acquiring boundary conditions and analyzing and calculating deformation of each element of the analysis model, and a mapping function for reflecting a state of each element in a planar shape to a product shape. The computer is used to calculate the linear displacement of each node in a direction perpendicular to the plane and the direction of the normal vector of each node when the analysis model created by the analysis model creation function is developed on a predetermined plane. Displacement calculation means for calculating the amount of rotational displacement required to make the rotation coincide with the direction orthogonal to the plane, and the analysis calculation by the analysis function described above. The method is characterized in that the linear displacement amount and the rotational displacement amount calculated by the displacement amount calculating unit are functioned as a displacement amount dividing unit that divides the linear displacement amount and the rotational displacement amount into a plurality of pieces in order to be a boundary condition of each step in the case of performing the steps separately. I do.
[0025]
According to a seventh aspect of the present invention, in the program according to the sixth aspect, when the computer is caused to function as step number setting means for setting the number of steps of analysis calculation and to function as displacement amount dividing means. The amount of linear displacement and the amount of rotational displacement are divided into the same number as the number of steps set by the function as the number-of-steps setting means.
[0026]
According to an eighth aspect of the present invention, in the program according to the sixth or seventh aspect, when the computer is caused to function as the displacement amount dividing means, a straight line of each node calculated by the displacement amount calculating means is used. The displacement amount and the rotational displacement amount are equally divided.
[0027]
According to a ninth aspect of the present invention, in the program according to the sixth or seventh aspect, the computer functions as a displacement characteristic setting means for setting a characteristic of each node with respect to a time course of a linear displacement and a rotational displacement. And, when functioning as displacement amount dividing means, dividing the linear displacement amount and rotational displacement amount of each node calculated by the displacement amount calculating means based on the displacement characteristics set by the displacement characteristic setting means. Features.
[0028]
The invention according to claim 6 to claim 9 corresponds to the invention relating to the system according to claim 1 to claim 4, wherein the program is stored in an analysis model creation function, an analysis function, and a mapping function. By applying the present invention to a computer provided with the above-described configuration, the same operation as the system according to the first to fourth aspects can be obtained.
[0029]
Further, the invention according to claim 10 relates to a computer-readable recording medium on which the program according to any one of claims 6 to 9 is recorded. By applying the present invention to a computer having a function, an analysis calculation function, a mapping function, and the like, the same operation as the system according to any one of claims 1 to 4 can be obtained.
[0030]
BEST MODE FOR CARRYING OUT THE INVENTION
Next, a press molding simulation system according to an embodiment of the present invention will be described.
[0031]
FIG. 1 shows the configuration of a computer serving as the center of this system. This computer 10 includes a central processing unit 11, an input device 12 used for setting various conditions and controlling the system, and a recording medium 20. A data reading device 13, a recording device 14 provided with a program recording unit 14a and a data file recording unit 14b, a display device 15 for displaying an input screen and calculation results, and a printing device 16 for printing calculation results and the like. Having.
[0032]
As shown in FIG. 2, a main program and a plurality of subprograms are recorded in the program recording unit 14a of the recording device 14, and the subprograms include an analysis model creation subprogram, a basic data creation subprogram, an analysis calculation A subprogram, a mapping processing subprogram, and a drawing subprogram are recorded.
[0033]
As shown in FIG. 3, an analysis model data file, an analysis basic data file, a material attribute data file, and an analysis result data file are recorded in the data file recording unit 14b of the recording device 14. It has become.
[0034]
Among these data files, the analysis model data file is composed of a table in which various data on the analysis model obtained by dividing the product to be analyzed into finite elements are recorded, and the coordinate data of each node shown in FIG. 4 is recorded. 5 includes a node coordinate table, an element configuration table in which nodes constituting each element shown in FIG. 5 are recorded, and a node configuration table in which elements constituting each node shown in FIG. 6 are recorded.
[0035]
Prior to the analysis calculation, the analysis basic data file is composed of a table in which various data calculated or set in advance are recorded. The element vector table shown in FIG. 7, the node vector table shown in FIG. , A rotational displacement table shown in FIG. 10, a constraint condition table shown in FIG. 11, a displacement characteristic curve table shown in FIG. 12, and the like.
[0036]
The material attribute data file is composed of a table shown in FIG. 13 which records various material attributes of the analysis target product necessary for the analysis calculation.
[0037]
Further, as shown in FIG. 14, the analysis result data file is configured by a table for recording analysis result data such as stress, strain, and thickness for all elements of the analysis model for each of a plurality of analysis steps.
[0038]
Next, the operation of the system will be described with reference to the flowchart of FIG. 15 showing the operation of the main program recorded in the program recording section 14a of the recording device 14.
[0039]
First, as a process P1, CAD data of a product to be analyzed is read from the recording medium 20 by the data reading device 13 of the computer 10 shown in FIG. At this time, a screen 15a displaying the entire shape of the product is displayed on the display device 15 of the computer 10, as shown in FIG.
[0040]
Next, as a process P2, an analysis model in which the shape of the product is divided into finite elements based on the CAD data is created by the analysis model creation subprogram recorded in the program recording unit 14a. At this time, as shown in FIG. 17, the created analysis model is displayed on the display device 15 as a screen 15b, and as shown in FIG. 18, the node numbers N1 and N2 are assigned to the nodes of the mesh forming the analysis model. , And element numbers E1, E2,.
[0041]
Then, the x-, y-, and z-coordinates of the nodes N1, N2,... When the analysis model is placed in a predetermined xyz coordinate system are read by the analysis model creation subprogram, and are stored in the node coordinate table of FIG. Written for each node. Further, the nodes constituting each of the elements E1, E2,..., That is, the numbers of the nodes around the element are read and written into the element configuration table shown in FIG. 5 for each element. Further, each of the nodes N1, N2. The constituent elements, that is, the numbers of the elements surrounding the node are read and written into the node configuration table of FIG. 6 for each node.
[0042]
Next, as process P3, analysis basic data for the analysis model is created by a basic data creation subprogram recorded in the program recording unit 14a.
[0043]
Here, the operation of this subprogram will be described with reference to the flowchart shown in FIG.
[0044]
First, in process P21, the coordinate data of the nodes, element configuration data, and node configuration data obtained by the above-described analysis model creation processing are read from the tables of FIGS. 4, 5, and 6, respectively. Next, in process P22, The normal vectors (hereinafter, referred to as “element vectors”) VE1, VE2,... For each of the elements E1, E2,.
[0045]
The method of calculating the element normal vectors VE1, VE2,... Will be described with respect to the element E1, for example. First, the fact that the element E1 is composed of nodes N1, N2, N5, N4 is read from the element configuration table of FIG. The x-coordinate, y-coordinate, and z-coordinate of the nodes N1, N2, N5, and N4 are read from the node coordinate table of FIG. 4, and the inclination of the element E1 when the element E1 is a plane is obtained based on these coordinate data. Then, the x component, the y component, and the z component of the unit vector in the direction perpendicular to the plane are obtained, and the vector having these components is set as the element vector VE1 of the element E1.
[0046]
In this way, the x, y, and z components of the element vectors VE1, VE2,... For each of the elements E1, E2,.
[0047]
Next, in process P23, normal vectors (hereinafter, referred to as “node vectors”) VN1, VN2,... For each of the nodes N1, N2,.
[0048]
The method of calculating the node normal vectors VN1, VN2,... Will be described with respect to the node N5 shown in FIG. 18, for example. Elements constituting the node N5, that is, elements surrounding the node N5 are read from the node configuration table of FIG. In this case, in the example of FIG. 18, the elements E1, E2, E4, and E3 are read.
[0049]
Then, the x, y, and z components of the element vectors VE1, VE2, VE4, and VE3 for these elements E1, E2, E4, and E3 are read from the element vector table of FIG. A unit vector having the average value of each component as the x component, the y component, and the z component is obtained, and this is set as a node vector VN5 for the node N5.
[0050]
In this manner, the node vectors VN1, VN2,... For each of the nodes N1, N2,... Are sequentially obtained, and the x-, y-, and z-components are written in the node vector table of FIG.
[0051]
Next, in process P24, from the nodal vectors VN1, VN2... For all the nodes N1, N2..., As the average direction of the directions of these vectors, the developing direction of the analytical model, Then, an XYZ coordinate system in which the direction is, for example, the Z direction is set as shown in FIG. Then, a plane orthogonal to the development direction Z is set on this XYZ coordinate system, and this is used as a virtual development plane P of the analysis model. In this case, the position of the development plane P in the Z direction can be set arbitrarily.
[0052]
Next, in process P25, linear displacement amounts D1, D2,... In the development direction Z of the nodes N1, N2,... When the analysis model is developed on the virtual development plane P are calculated, and are calculated as straight lines shown in FIG. Write in the displacement amount table together with the development direction.
[0053]
In this case, the calculation of the linear displacement amounts D1, D2,... Is performed by using the coordinates in the xyz coordinate system of the nodes N1, N2,. After conversion into coordinates in the coordinate system, the XYZ coordinate system is used to determine the distance in the Z direction from each of the nodes N1, N2... To the virtual development plane P. FIG. The quantity D1 is illustrated.
[0054]
In the process P26, the rotational displacement amounts (AX1, AY1, AZ1) around the axes X, Y, Z of the node vectors VN1, VN2... When the analysis model is developed on the virtual development plane P, (AX2) , AY2, AZ2)... Are written in the rotational displacement amount table of FIG.
[0055]
This rotational displacement amount is a rotation angle around each axis required to make the direction of the node vector coincide with the development direction of the analysis model. FIG. 20 illustrates the node vector VN1 of the node N1 as an example. The rotation angle around the Y axis required to move VN1 on a plane parallel to the YZ plane is AY1, and the rotation angle around the X axis required to further move this vector on the plane parallel to the ZX plane. Is AX1, the direction of the node vector VN1 coincides with the Z direction due to these rotation angles.
[0056]
That is, these rotation angles AY1 and AX1 are the amounts of rotational displacement around the Y axis and the X axis for making the vector VN1 coincide with the development direction. In this case, since the developing direction and the Z direction match, the amount of rotational displacement around the Z axis is zero.
[0057]
Next, in process P27, for each of the nodes N1, N2,..., A displacement characteristic curve of a linear displacement and a rotational displacement is defined, and data indicating the curve is recorded in the displacement characteristic curve table of FIG.
[0058]
This displacement characteristic curve shows the characteristic of the change over time of the linear displacement and the rotational displacement when the analysis model is developed on the virtual development plane P, and is set as follows, for example.
[0059]
As shown in FIG. 21, it is assumed that, for example, the node Na of the analysis model is located at a position indicated by reference numeral a with respect to the node Nb, and is displaced to a position indicated by reference numeral a when it is developed into a planar shape. At this time, the final linear displacement amount of the node Na is D, and the rotational displacement amount of the node vector VNa around an axis orthogonal to the paper surface is A. Further, assuming that the linear displacement amount of the node Na up to the middle stage of the deformation process indicated by the symbol c is d, the rotational displacement amount of the node vector VNa is a, and the distance between the nodes Na and Nb is R, the deformation process is The relationship between the linear displacement d and the rotational displacement a of
d = DR-sin (A-a)
R = D / sinA
Because
d = D · [1-sin (A−a) / sinA]
It becomes.
[0060]
Here, for example, the final linear displacement amount D is set to 3, the final rotational displacement amount A is set to 0.785 (45 °), and a value obtained by dividing the final rotational displacement amount A into 10 as the intermediate stage rotational displacement amount a, By setting 0.079, 0.157, and the like, and calculating the linear displacement amount at each stage for these values, the values shown in the table of FIG. 22 are obtained.
[0061]
Further, for example, assuming that the rotational displacement changes linearly with the passage of time, if the above calculation result is shown as a characteristic with respect to the passage of time, a linear displacement characteristic curve LD and a rotational displacement characteristic curve LA as shown in FIG. Is obtained.
[0062]
Here, in practice, it is necessary to define the characteristics of the linear displacement and the rotational displacement around the X-axis and the characteristic of the rotational displacement around the Y-axis with respect to the passage of time. It is also possible to calculate and define what characteristics the rotational displacement around the X-axis and the Y-axis at that time changes.
[0063]
In addition, since there are many nodes that are deformed with the same characteristics over time even if the final displacement amount is different, when defining the displacement characteristics, the final displacement amount is set to 1 and an intermediate stage for each time is set. Is defined as a ratio to the total displacement, and when this is used as the displacement characteristic of a specific node, the value (ratio) on the characteristic curve is multiplied by the final displacement of that node to The intermediate displacement may be obtained. In this way, the displacement characteristic curve is shared for a plurality of nodes having the same displacement characteristic, so that the setting work and the recording capacity can be reduced.
[0064]
Specifically, as shown in the displacement characteristic curve table of FIG. 11, using the characteristic curve number as a key, the linear displacement up to each time when the analysis time is divided into ten equal parts, the rotational displacement around the X axis, and the rotation around the Y axis And a linear displacement amount table in which the linear displacement amount and the rotational displacement amount at each of the nodes in FIGS. 9 and 10 are recorded, and the rotational displacement amount is recorded as a ratio when the final displacement amount is set to 1. In the table, record the number of the displacement characteristic curve applied to each node.
[0065]
Next, in processes P28 and P29, two nodes for which displacements in the X direction and the Y direction are prohibited when the analysis model is developed are specified as constraint conditions at the time of analysis calculation. This is because, when the analysis model is developed into a planar shape, if all the nodes are movable in the X and Y directions, the position of the entire model on the XY plane is not determined and calculation becomes impossible. This is to prevent it.
[0066]
For example, if the node Nx is designated as a node whose displacement in the X direction is prohibited and the node Ny is designated as a node whose displacement in the Y direction is prohibited, the node numbers Nx, It is recorded that Ny and the displacements of the node in the X direction and the Y direction are each 0.
[0067]
Thus, the basic data creation processing by the basic data creation subprogram ends, and in process P30, the obtained data is stored in the analysis basic data file. Then, the operation of the computer 10 returns to the processing by the main program shown in the flowchart in FIG. 15, and executes the process P4 and thereafter.
[0068]
In the process P4 of the main program, the number of analysis calculation steps is set based on the input from the input device 12 of the computer 10. In the process P5, data relating to the material of the product to be analyzed this time, that is, data necessary for analyzing the Young's modulus, Poisson's ratio, yield stress, tangent coefficient, and the like of the material are read from the table of FIG. Then, the analysis processes of the processes P6 to S10 are started.
[0069]
In this analysis process, first, in process P6, boundary conditions necessary for the analysis calculation in the first step are created.
[0070]
That is, from the tables in FIGS. 9 and 10, the final linear displacements D1, D2,... Of the nodes N1, N2... When the analysis model is expanded to the planar shape, and the X-axis of the node vectors VN1, VN2. , AY1, AY2,..., And the displacement characteristic curve numbers LD1,..., LAX1,.
[0071]
Assuming that the number of analysis steps set in the process P5 is, for example, 10, the analysis time of each displacement characteristic curve defined using the numbers LD1,..., LAX1,. The values (ratio) at the first time when the data are divided into ten equal parts are read out, and the read values and the linear displacement amounts D1, D2,..., And the rotational displacement amounts AX1, AX2,. Is multiplied for each node to calculate the boundary condition in the first step of each node.
[0072]
Thereby, the boundary conditions in the first step of the linear displacement of each of the nodes N1 and N2 and the boundary conditions in the first step of the rotational displacement of each of the node vectors VN1, VN2... Around the X axis and the Y axis are obtained. become.
[0073]
Here, in the above example, since the characteristic curve of each displacement is defined by ten pieces of data and the set number of analysis steps is also ten, the data used in the definition of the characteristic curve is directly used to calculate the boundary condition. However, if they do not match, the boundary condition of each step is calculated from the characteristic curve using an interpolation method.
[0074]
When the boundary conditions of the linear displacement, the rotational displacement around the X-axis, and the rotational displacement around the Y-axis for all the nodes are obtained in this way, then, in the process P7 of the main program, each of the nodes N1, N1 shown in FIG. The initial conditions are set based on the values obtained by converting the x, y, and z coordinates of the N2... On the XYZ coordinate system, and the boundary conditions obtained as described above and the boundary conditions of FIG. Using the constraint conditions for the two nodes Nx and Ny shown in the table, the analysis calculation of the first step is executed by the analysis calculation subprogram.
[0075]
This calculation is performed by a static implicit method, a rigidity matrix is set using material attribute data such as Young's modulus and Poisson's ratio, and the stress, strain, and thickness of each element E1, E2... Calculate changes such as In the process P8, the calculation result is recorded in the table of FIG. 14 as an analysis result for all the elements E1, E2,.
[0076]
Next, in a process P9, it is determined whether or not the analysis calculation of all the steps has been completed. At this time, since the analysis calculation has not been completed yet, the process P10 is executed, and the first process recorded in the analysis result table in the process P8 is performed. The analysis result of the step is read out, and based on the result, the initial condition of the analysis calculation of the second step is set, and the process P6 is executed again. Create boundary conditions.
[0077]
In this case, the final linear displacement amounts D1, D2,... Of the nodes N1, N2,... Read from the tables of FIG. , AY1, AY2,... Are multiplied by the values (ratio) at the second time when the analysis time of each displacement characteristic curve defined in the table of FIG. Thus, the boundary condition in the second step of each node or the node vector is obtained.
[0078]
Then, using these boundary conditions, the initial conditions set based on the previous analysis result set in the process P10, and the constraint conditions recorded in the table of FIG. 11, the analysis calculation subprogram executes the process P7. The analysis calculation is executed again, and in the process P8, the calculation result is recorded in the analysis result table of FIG. 14 as the analysis result for all the elements E1, E2,.
[0079]
Similarly, the processes P6 to P10 are repeated until it is determined that the analysis calculation of all the steps is completed in the process P9, and the initial condition set based on the previous analysis result and the analysis time of the displacement characteristic curve are set to the number of steps. The analysis calculation is repeatedly executed using the boundary condition calculated using the value at each time divided by the above and the constraint condition recorded in the table, and the obtained calculation result is calculated for each step in FIG. Record in the analysis result table.
[0080]
When the analysis calculation for the set number of steps is completed in this way, the analysis model has been developed into a planar shape, and then the processes P9 to P11 are executed to perform the mapping process.
[0081]
This processing is executed by the mapping processing subprogram, and the calculation result of the final step recorded in the analysis result table of FIG. 14, that is, the stress, strain, and thickness of each element E1, E2. By applying data such as height to the elements E1, E2,... Of the analysis model before development, distribution of stress, strain, thickness, etc. of each part when a uniform plate-shaped blank is pressed into the product shape. Is what you want.
[0082]
Then, in the process P12, the distribution of the thickness on the product shape obtained by the mapping process is displayed by the drawing subprogram, for example, as shown on a screen 15c in FIG.
[0083]
As described above, the state of the press forming from the sheet material to the product shape given by the CAD data is simulated, and it is possible to infer in advance defects such as wrinkles and cracks during the press forming. .
[0084]
In the above embodiment, the characteristics of the linear displacement of each node and the rotational displacement around the X axis and the Y axis with respect to the passage of time are defined, and this is divided to obtain the boundary condition of each step. However, if simplification of the system configuration or quick calculation is required, both the linear displacement and the rotational displacement may be processed as if they change linearly with time. In this case, the boundary condition in each step is a value when the final displacement amount is equally divided by the set number of steps.
[0085]
Further, as described above, the system according to this embodiment is configured to operate by a plurality of programs, among which a finite element division program, an analysis calculation program, a mapping processing program, a drawing program, etc. Existing ones can be used. Therefore, simply loading the basic data creation program and the program for associating these programs into a general-purpose system in which these existing programs are installed makes it possible to operate the system as a system for implementing the present invention. This makes it possible to create basic data for each product in advance using a basic data creation program.
[0086]
In this case, the newly read basic data creation program and the like constitute an embodiment of the simulation program according to the present invention, and the computer-readable recording medium on which the program is recorded is the recording medium according to the present invention. This constitutes an embodiment of a medium.
[0087]
【The invention's effect】
As described above, according to the system, the program, and the recording medium according to the present invention, in addition to the calculation speed, which is a feature of the conventional one-step method, as a press molding simulation system using the finite element method, The accuracy of the simulation is improved. Therefore, a system that can be used in a wide range of situations from a product design stage to a mold manufacturing stage is realized.
[0088]
Further, since an existing analysis model creation program, analysis calculation program, and the like for dividing a product shape into finite elements can be used, there is an advantage that the above-described effect can be realized using an existing general-purpose system or the like. .
[Brief description of the drawings]
FIG. 1 is a block diagram showing an overall configuration of a system according to an embodiment of the present invention.
FIG. 2 is an explanatory diagram of a program recorded in a recording device of the system.
FIG. 3 is an explanatory diagram of a data file.
FIG. 4 is an explanatory diagram of a node coordinate table included in a finite element data file.
FIG. 5 is an explanatory diagram of an element configuration table.
FIG. 6 is an explanatory diagram of a node configuration table.
FIG. 7 is an explanatory diagram of an element vector table included in an analysis basic data file.
FIG. 8 is an explanatory diagram of a node vector table.
FIG. 9 is an explanatory diagram of a linear displacement amount table.
FIG. 10 is an explanatory view of a rotational displacement amount table.
FIG. 11 is an explanatory diagram of a displacement characteristic curve table.
FIG. 12 is an explanatory diagram of a constraint condition table.
FIG. 13 is an explanatory diagram of a table constituting a material attribute data file.
FIG. 14 is an explanatory diagram of a table constituting an analysis result data file.
FIG. 15 is a flowchart showing an operation of a main program for controlling the entire system.
FIG. 16 is an explanatory diagram of a screen for displaying a CAD drawing.
FIG. 17 is an explanatory diagram of a screen displaying an analysis model obtained by dividing a finite element.
FIG. 18 is an explanatory diagram of elements, nodes, and the like of an analysis model.
FIG. 19 is a flowchart showing the operation of an analysis basic data creation program.
FIG. 20 is an explanatory diagram of a displacement amount when an analysis model is developed.
FIG. 21 is an explanatory diagram illustrating an example of a method of defining a displacement characteristic curve.
FIG. 22 is a table showing displacement characteristic curve data obtained by the method.
FIG. 23 is a graph illustrating the displacement characteristic curve.
FIG. 24 is an explanatory diagram of a screen showing a simulation result.
[Explanation of symbols]
10 computer
11 Central processing unit
12 Input device
13 Data reading device
14 Storage device
15 Display device
