JP3897477B2 - Stress-strain relationship simulation method and springback amount prediction method - Google Patents

Description

ただし、
Ｓ_i （ｉ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に作用する垂直応力σ_i に関する偏差応力
Ｓ_ij（ｉ，ｊ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に、指標ｊに対応する座標軸の方向に作用するせん断応力σ_ijに関する偏差応力
α_i （ｉ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に作用する背応力
α_ij（ｉ，ｊ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に、指標ｊに対応する座標軸の方向に作用する背応力
Ｙ₀ ：材料の塑性変形が初めて開始されるときの降伏応力である初期降伏応力
Ｒ：材料の等方硬化量
なる式で表された(1) 項に記載の応力−ひずみ関係シミュレート方法。
(3) 前記背応力αを表すベクトル｛α｝が、
【数９】

ただし、
｛α₁ ｝：背応力のうち移動硬化の非線形性に依存する成分である第１背応力成分を表すベクトル
｛α₂ ｝：背応力のうち移動硬化の線形性に依存する成分である第２背応力成分を表すベクトル
なる式で表され、その第１背応力成分α₁ の増分ｄα₁ を表すベクトル｛ｄα₁ ｝が、
【数１０】

ただし、
Ｃ：非線形移動硬化の収束の速さを表す係数
２ａ／３：非線形移動硬化の収束値
｛ｄε^P ｝：塑性ひずみ増分ｄε^P _i （ｉ＝ｘ，ｙ，ｚ）を表すベクトル
ε^P _eq：相当塑性ひずみ
なる式で表され、前記第２背応力成分α₂ の増分ｄα₂ を表すベクトル｛ｄα₂ ｝が、
【数１１】

ただし、
Ｈ：線形移動硬化の大きさを表す係数
なる式で表された(2) 項に記載の応力−ひずみ関係シミュレート方法。
(4) さらに、前記等方硬化量と前記塑性ひずみとの関係を、塑性ひずみの大小を問わず等方硬化量が塑性ひずみに応じて増加するものとするとともに、それら等方硬化量と塑性ひずみとの関係を、取得された複数の実験値に基づいて同定する関係同定工程を含む(1) ないし(3) 項のいずれかに記載の応力−ひずみ関係シミュレート方法〔本項の一部が請求項３に相当する〕。
この方法によれば、塑性ひずみの大小を問わず等方硬化量が塑性ひずみに応じて増加するため、高ひずみ領域においても等方硬化量が塑性ひずみによって変化することとなり、バウシンガ効果を精度よくシミュレートするのに好都合である。
(5) 前記等方硬化量Ｒの増分ｄＲが、
【数１２】

ただし、
ｂ：等方硬化の大きさを表す定数
なる式で表された(4) 項に記載の応力−ひずみ関係シミュレート方法〔本項の一部が請求項４に相当する〕。
(6) 前記係数Ｃが、前記塑性ひずみに応じて減少する傾向を有する関数で定義された(3) ないし(5) 項のいずれかに記載の応力−ひずみ関係シミュレート方法〔本項の (3) 項に従属する部分から (2) 項および (3) 項の具体的な式の限定を除いたものが請求項１の一部に相当する〕。
(7) 前記係数Ｃを表す関数が、前記塑性ひずみとして相当塑性ひずみε^P _eqを用いた場合に、その相当塑性ひずみε^P _eqを１／√（ε^P _eq）という形態で組み込んだものである(6) 項に記載の応力−ひずみ関係シミュレート方法。
(8) 前記係数Ｃが、
【数１３】

ただし、
ε^P _eq：相当塑性ひずみ
Ｃ’，ε₀ ：定数
なる式で表された(3) ないし(7) 項のいずれかに記載の応力−ひずみ関係シミュレート方法〔請求項２の一部に相当する〕。
(9) 前記係数Ｈが、前記塑性ひずみに応じて増加する傾向を有する関数で定義された(3) ないし(8) 項のいずれかに記載の応力−ひずみ関係シミュレート方法〔本項の (3) 項に従属する部分から (2) 項および (3) 項の具体的な式の限定を除いたものが請求項１の一部に相当する〕。
(10)前記係数Ｈを表す関数が、前記塑性ひずみとして相当塑性ひずみε^P _eqを用いた場合に、その相当塑性ひずみε^P _eqを自然対数ｌｎ（ε^P _eq）または常用対数ｌｏｇ（ε^P _eq）という形態で組み込んだものである(9) 項に記載の応力−ひずみ関係シミュレート方法。
(11)前記係数Ｈが、
【数１４】

ただし、
ε₁ ：境界値
ε₂ ：境界値（＞ε₁ ）
Ｈ₁ ，ｐ，ｑ，Ｈ₂ ：定数
なる式で表された(3) ないし(10)項のいずれかに記載の応力−ひずみ関係シミュレート方法〔請求項２の一部に相当する〕。
(12)(1) ないし(11)項のいずれかに記載の応力−ひずみ関係シミュレート方法を実施するためにコンピュータにより実行されるプログラムをコンピュータ読み取り可能に記録した記録媒体。
ここに、「記録媒体」には例えば、フレキシブルディスク，磁気テープ，磁気ディスク，磁気ドラム，磁気カード，光ディスク，光磁気ディスク，ＲＯＭ，ＣＤ−ＲＯＭ，ＩＣカード，穿孔テープ等がある。この説明は後述の「記録媒体」についても適用される。
(13)(1) ないし(11)項のいずれかに記載の降伏関数を用いることにより、材料が塑性加工の終了後に弾性回復する量であるスプリングバック量を予測する方法であって、
前記材料と材質が実質的に同じ試験片を引張方向と圧縮方向との一方向に負荷して塑性変形させた後に除荷し、さらに、逆方向に負荷して塑性変形させるとともに、その間にその試験片の応力−ひずみ関係の実験値を複数の離散値として取得する実験値取得工程と、
前記降伏関数における前記背応力の前記係数を、前記塑性ひずみを変数とする係数関数により定義された変数とするとともに、その係数関数を、取得された複数の実験値に基づいて同定する係数関数同定工程と、
前記塑性加工の終了時に前記材料に残存する応力−ひずみ関係を予測する残存応力−ひずみ関係予測工程と、
その予測された残存応力−ひずみ関係と、同定された係数関数により定義された前記係数を含む前記降伏関数とに基づき、前記スプリングバック量を予測するスプリングバック量予測工程と
を含むことを特徴とするスプリングバック量予測方法〔本項の一部が請求項５に相当する〕。
プレス成形等、塑性加工においては、材料が負荷された後に除荷されるが、除荷後に材料にスプリングバックが生じ、材料の除荷終了時における形状（最終形状）が、負荷終了時（除荷開始時）における形状とは異なったものとなる。そのため、材料の塑性加工後の形状精度を確保するために、プレス型等、塑性加工のための工具の形状を、材料の目標形状に対応する工具形状にスプリングバック量を見込んだ工具形状に修正することが必要となる。この工具修正は、形状が暫定的である工具を用いて塑性加工を試行して材料のスプリングバック量を実測することによって取得し、その実測値を考慮して工具形状を修正することによって実現可能であるが、塑性加工の試行を行うことなくスプリングバック量を高い精度で予測することができれば、工具修正に必要な時間が短縮される。工具修正に必要な時間が短縮されれば、生産準備リードタイムが短縮され、生産効率が向上する。しかし、従来、材料のスプリングバック量を精度よく予測することが困難であった。前述のように、従来の手法では、高ひずみ領域において材料のバウシンガ効果を精度よくシミュレートすることができなかったからである。
これに対して、本項に記載のスプリングバック量予測方法によれば、高ひずみ領域においても材料のバウシンガ効果を精度よくシミュレートすることが可能となる。一方、材料のバウシンガ効果のシミュレート精度が高くなれば、スプリングバック量の予測精度も高くなる。したがって、この方法によれば、スプリングバック量を精度よく予測することが可能となる。
(14)さらに、前記等方硬化量と前記塑性ひずみとの関係を、塑性ひずみの大小を問わず等方硬化量が塑性ひずみに応じて増加するものとするとともに、それら等方硬化量と塑性ひずみとの関係を、取得された複数の実験値に基づいて同定する関係同定工程を含み、かつ、前記スプリングバック量予測工程が、(a) 予測された残存応力−ひずみ関係と、(b) 同定された係数関数により定義された前記係数と、同定された関係を表す係数とを含む前記降伏関数とに基づき、前記スプリングバック量を予測するものとされた(13)項に記載の応力−ひずみ関係シミュレート方法。
この方法によれば、塑性ひずみの大小を問わず等方硬化量が塑性ひずみに応じて増加するため、高ひずみ領域においても等方硬化量が塑性ひずみによって変化することとなり、バウシンガ効果を精度よくシミュレートするのに好都合である。
(15)(13)または(14)項に記載のスプリングバック量予測方法を実施するためにコンピュータにより実行されるプログラムをコンピュータ読み取り可能に記録した記録媒体〔本項の一部が請求項６に相当する〕。
【０００７】
【発明の実施の形態】
以下、本発明のさらに具体的な実施の形態の一つを図面に基づいて詳細に説明する。
【０００８】
本発明の一実施形態は、応力−ひずみ関係シミュレート方法を含むスプリングバック量予測方法であり、図１には、そのスプリングバック量予測方法を実施するのに好適なスプリングバック量予測装置（以下、単に「予測装置」という）が示されている。
【０００９】
予測装置は、入力装置１０とコンピュータ１２と出力装置１４と外部記憶装置１６とを備えている。入力装置１０はマウス，キーボード等を含むように構成される。コンピュータ１２はＣＰＵ等、プロセッサ２０と、ＲＯＭ，ＲＡＭ，ハードディスク等、メモリ２２と、それらプロセッサ２０とメモリ２２とを接続するバス２４とを含むように構成される。出力装置１４はディスプレイ，プリンタ，プロッタ等を含むように構成される。外部記憶装置１６は、ＣＤ−ＲＯＭ３０，書き込み可能なＦＤ３２（フレキシブルディスク）等、記録媒体が装填可能となっていて、装填状態においては、記録媒体に対するデータの読み取りおよび書き込みが必要に応じて行われる。
【００１０】
本実施形態においては、上記スプリングバック量の予測に必要なプログラムがＣＤ−ＲＯＭ３０に記憶され、また、その予測に必要なデータがＦＤ３２に記憶されている。スプリングバック量の予測に必要なプログラムには、塑性構成式同定プログラムと、成形シミュレーションプログラムと、スプリングバック量予測プログラムとがある。そして、スプリングバック量の予測時には、それらＣＤ−ＲＯＭ３０およびＦＤ３２から必要なプログラムおよびデータが読み出されてコンピュータ１２のＲＡＭまたはハードディスクに転送され、その後、プロセッサ２０によりそのプログラムが実行される。
【００１１】
図２には、本実施形態であるスプリングバック量予測方法が工程図で表されている。
【００１２】
概略的に説明すれば、スプリングバック量予測方法は、鋼板である材料をプレス成形した後にその材料に生じるスプリングバック量を予測するものである。そのプレス成形は、塑性ひずみが０．１５，０．２０というように、比較的高い塑性ひずみを有する部分が生じるように行われる。また、このスプリングバック量予測方法に含まれる応力−ひずみ関係シミュレート方法は、プレス成形されるべき材料と材質が実質的に同じ試験片を用いることにより、その試験片の応力−ひずみ関係の実験値を複数の離散値として取得し、その取得した複数の実験値に基づき、その試験片の応力−ひずみ関係を構成する方程式である塑性構成式を同定するものである。この塑性構成式は、応力−ひずみ関係を連続的に記述するものである。また、スプリングバック量予測方法は、有限要素法を用いることにより、プレス成形直後に材料に残存する応力−ひずみ関係を予測するために、材料のプレス成形をシミュレートする成形シミュレーションを行い、その後、有限要素法を用いるとともに、同定された塑性構成式を用いることにより、材料のスプリングバック量を予測するものである。
【００１３】
さらに具体的に説明すれば、本実施形態においては、同図に示すように、まず、ステップＳ１において、作業者により、プレス成形されるべき材料と材質が実質的に同じ試験片が準備され、その後、試験機を用いることにより、その試験片に対して一軸応力状態において引張−圧縮試験が行われる。この試験は具体的には、試験片を引張方向に負荷して塑性変形させた後に除荷し、さらに、逆方向すなわち圧縮方向に負荷して塑性変形させるというものであり、その間にその試験片の応力−ひずみ関係の実験値が複数の離散値として取得される。
【００１４】
次に、ステップＳ２において、コンピュータ１２により、上記塑性構成式同定プログラムが実行され、それにより、上記取得された複数の実験値に基づいて上記塑性構成式が同定される。この同定については後に詳述する。
【００１５】
その後、ステップＳ３において、作業者により、上記成形シミュレーションを実行するのに必要なデータ、すなわち、数値解析用データが作成される。数値解析用データには、図３に示すように、材料を複数のメッシュに分割するためのメッシュデータと、プレス成形時に材料に付与される応力に関する条件と、複数のメッシュに関する境界条件とを含んでいる。
【００１６】
続いて、ステップＳ４において、コンピュータ１２により、上記成形シミュレーションプログラムが実行され、それにより、プレス成形直後に材料に残存する応力−ひずみ関係が予測される。さらに、プレス成形直後の材料の形状も予測される。この成形シミュレーションプログラムの一例は、ＬＳＴＣ社により製造され、「ＬＳ−ＤＹＮＡ３Ｄ」という名称で日本総合研究所により販売されたものである。
【００１７】
その後、ステップＳ５において、上記ステップＳ２において同定された塑性構成式により、上記スプリングバック量予測プログラムにおいて、スプリングバック量を予測すべき材料の応力−ひずみ関係が定義される。なお、このスプリングバック量予測プログラムは、材料の応力−ひずみ関係を任意に定義可能に設計されている。
【００１８】
続いて、ステップＳ６において、コンピュータ１２により、上記スプリングバック量予測プログラムが実行され、それにより、スプリングバック量が予測される。このスプリングバック量予測プログラムの一例は、日本総合研究所により製造され、同研究所により「ＪＯＨ−ＮＩＫＥ３Ｄ」という名称で販売されたものである。
【００１９】
以上で、スプリングバック量予測方法の一回の実行が終了する。なお、図３には、それら３つのプログラム相互の関係が概念的に示されている。
【００２０】
ここで、塑性構成式の同定について詳しく説明する。
【００２１】
１．従来の複合硬化モデル
雑誌「Journal of Pressure Vessel Technology (MAY 1983, Vol. 105/154-159)」には、J.L.ChabocheとG.Rousselierとにより著された、「On the Plastic and Viscoplastic Constitutive Equations」と題する論文が掲載されている。この論文においては、弾塑性体の応力−ひずみ関係を与える構成方程式である塑性構成式が記載されている。この塑性構成式は、材料の応力−ひずみ関係を、等方硬化モデルと移動硬化モデルとを組み合わせた複合硬化モデルで近似することにより、材料のバウシンガ効果を表現し得るものである。
【００２２】
ここに「等方硬化」とは、加工硬化により降伏曲面がその中心位置およびその形状が変化することなく等方的に拡大することを意味する。また「移動硬化」とは、加工硬化により降伏曲面がその形状および大きさが変化することなくその中心位置が移動することを意味する。また「降伏曲面」とは、降伏関数を応力空間または応力平面に表した曲面（多角形を含む）をいい、「降伏関数」とは、材料が降伏するときの応力成分σ_ij（ｉ，ｊ＝ｘ，ｙ，ｚ）の関数をいう。
【００２３】
Chaboche-Rousselier により提案された塑性構成式としての降伏関数ｆは、次の式(1) で表される。
【００２４】
【数１５】

【００２５】
ただし、
Ｓ_i （ｉ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に作用する垂直応力σ_i に関する偏差応力
Ｓ_ij（ｉ，ｊ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に、指標ｊに対応する座標軸の方向に作用するせん断応力Ｓ_ij（＝τ_ij）に関する偏差応力α_i （ｉ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に作用する背応力
α_ij（ｉ，ｊ＝ｘ，ｙ，ｚ）：指標ｉに対応する座標軸に垂直な面に、指標ｊに対応する座標軸の方向に作用する背応力
Ｙ₀ ：材料の塑性変形が初めて開始されるときの降伏応力である初期降伏応力
Ｒ：材料の等方硬化量
【００２６】
なお、上記式(1) は、ミーゼス型の降伏条件と等方硬化と線形硬化とを前提として成立している。
【００２７】
また、この式(1) において、背応力α_i ，α_ijは、降伏曲面の中心の移動量を表すことから、移動硬化モデルにおける移動硬化量ということができる。なお、図４には、引張−圧縮時における応力−ひずみ関係がグラフで表されるとともに、背応力αと等方硬化量Ｒと初期降伏応力Ｙ₀ との関係が示されている。すなわち、この図には、複合硬化モデルが概念的にグラフで示されているのである。
【００２８】
一軸応力状態においては、Chaboche-Rousselier により、背応力α_i （ｉ＝ｘ，ｙ，ｚ）を表すベクトル｛α｝は、次の式(2) により定義される。
【００２９】
【数１６】

【００３０】
ただし、
｛α₁ ｝：背応力のうち移動硬化の非線形性に依存する成分である第１背応力成分を表すベクトル
｛α₂ ｝：背応力のうち移動硬化の線形性に依存する成分である第２背応力成分を表すベクトル
【００３１】
ここに、第１背応力成分α₁ を表すベクトル｛α₁ ｝は、次の式(3) により定義される。
【００３２】
【数１７】

【００３３】
ただし、
Ｃ：非線形移動硬化の収束の速さを表す係数（定数）
２ａ／３：非線形移動硬化の収束値（定数）
｛ｄε^P ｝：塑性ひずみ増分ｄε^P _i （ｉ＝ｘ，ｙ，ｚ）を表すベクトル
ｄε^P _eq：相当塑性ひずみ増分
【００３４】
また、第２背応力成分α₂ を表すベクトル｛α₂ ｝は、次の式(4) により定義される。
【００３５】
【数１８】

【００３６】
ただし、
Ｈ：線形移動硬化の大きさを表す係数（定数）
【００３７】
また、等方硬化量Ｒの増分ｄＲは、次の式(5) により定義される。
【００３８】
【数１９】

【００３９】
ただし、
ｂ：等方硬化の大きさを表す係数（定数）
Ｑ：等方硬化の収束値（定数）
【００４０】
２．従来の複合硬化モデルの問題点
図５には、塑性ひずみε^P が０．０２より小さい低ひずみ領域に移行させるべく、材料を引っ張って塑性変形させた後に除荷し、さらにその材料を圧縮して塑性変形させた場合の応力−ひずみ関係が示されている。応力−ひずみ関係の実験値は複数の白丸が並んだグラフで表され、一方、Chaboche-Rousselier により提案された従来の複合硬化モデルに基づく予測値は破線のグラフで表されている。それらグラフから明らかなように、低ひずみ領域においては、従来の複合硬化モデルに基づく予測値は実験値と精度よく一致しており、バウシンガ効果が精度よくシミュレートできている。
【００４１】
図６には、塑性ひずみが０．０２より大きく、０．１５４より小さい高ひずみ領域に移行させるべく、材料を引っ張って塑性変形させた後に除荷し、さらにその材料を圧縮して塑性変形させた場合の応力−ひずみ関係が示されている。実験値は複数の白丸が並んだグラフで表され、一方、Chaboche-Rousselier により従来の複合硬化モデルに基づく予測値は破線のグラフで表されている。それらグラフから明らかなように、高ひずみ領域においては、従来の複合硬化モデルに基づく予測値は実験値と精度よく一致しておらず、バウシンガ効果が精度よくシミュレートできていない。
【００４２】
３．本実施形態における複合硬化モデル
(1) モデルの設計目標
低ひずみ領域のみならず高ひずみ領域においてもバウシンガ効果が精度よくシミュレートできるようにする。
【００４３】
(2) 複合硬化モデルの概要
▲１▼ 係数Ｃの変数化
前記式(3) における係数Ｃを、定数ではなく、相当塑性ひずみε^P _eqを変数とする関数（係数関数の一例である）で定義された変数とする。係数Ｃは、相当塑性ひずみε^P _eqに応じて減少する傾向を有する。係数Ｃは具体的には、次の式(6) で表される。
【００４４】
【数２０】

【００４５】
ただし、
Ｃ’：応力−ひずみ関係の実験値を用いて同定された係数（後述する）
ε₀ ：相当塑性ひずみε^P _eqが０である場合に上記式(6) の分母が０になってしまうことを防止するための定数（例えば、１０^-6程度の値）
【００４６】
▲２▼ 係数Ｈの変数化
前記式(4) における係数Ｈを、定数ではなく、相当塑性ひずみε^P _eqを変数とする関数（係数関数の一例である）で定義された変数とする。係数Ｈは、相当塑性ひずみε^P _eqに応じて増加する傾向を有する。係数Ｈは具体的には、次の式(7) ないし(9) で表される。
【００４７】
【数２１】

【００４８】
ただし、
ε₁ ：境界値
Ｈ₁ ｐ，ｑ，Ｈ₂ ：応力−ひずみ関係の実験値を用いて同定された定数
ε₂ ：境界値（＞ε₁ ）
【００４９】
▲３▼ 等方硬化量Ｒの構成式の変形
相当塑性ひずみε^P _eqの大小を問わず等方硬化量Ｒが相当塑性ひずみε^P _eqに応じて増加するようにし、それにより、低ひずみ領域のみならず高ひずみ領域においても等方硬化量Ｒが相当塑性ひずみε^P _eqに応じて増加するようにする。このようにすることは、高ひずみ領域においてもバウシンガ効果を精度よくシミュレートするのに好都合である。具体的には、前記式(5) を、次の式(10)に変形する。
【００５０】
【数２２】

【００５１】
その結果、等方硬化量Ｒと相当塑性ひずみε^P _eqとの間に、その相当塑性ひずみε^P _eqの大小を問わず、比例関係が成立することになる。
【００５２】
▲５▼ 係数Ｃ’の決定
前記式(6) における係数Ｃ’は、定数ではなく、引張終了時（除荷開始時）における相当塑性ひずみε^P _eqである基準ひずみε^P _T を変数とする関数で定義された変数とする。よって、係数Ｃ’は、各回の引張−圧縮状態ごとに一義的に決まってしまい、同じ回の引張−圧縮状態においては、相当塑性ひずみε^P _eqが変化しても変化しない。
【００５３】
係数Ｃ’を表す関数は例えば、次の式(11)で表される。
【００５４】
【数２３】

【００５５】
ただし、
Ｃ’₁ ，Ｃ’₀ ：応力−ひずみ関係の実験値を用いて同定された定数
【００５６】
なお、係数Ｃ’は、引張時には、基準ひずみε^P _T が０であるため、係数Ｃ’₀ と等しくなる。
【００５７】
また、係数Ｃ’は、他の式により定義可能であり、例えば、ｎ（≧２）次の基準ひずみε^P _T の項を含む多項式により定義可能である。
【００５８】
▲６▼ 係数ａの決定
前記式(3) における係数ａを、定数ではなく、基準ひずみε^P _T を変数とする関数で定義された変数とする。よって、係数ａは、上記係数Ｃ’と同様、各回の引張−圧縮状態ごとに一義的に決まってしまい、同じ回の引張−圧縮状態においては、相当塑性ひずみε^P _eqが変化しても変化しない。
【００５９】
係数ａを表す関数は例えば、次の式(12)で表される。
【００６０】
【数２４】

【００６１】
ただし、
ａ₁ ，ａ₀ ：応力−ひずみ関係の実験値を用いて同定された定数
【００６２】
なお、係数ａは、引張時には、基準ひずみε^P _T が０であるため、係数ａ₀ と等しくなる。
【００６３】
また、係数ａは、他の式により定義可能であり、例えば、ｎ（≧２）次の基準ひずみε^P _T の項を含む多項式により定義可能である。
【００６４】
(3) 引張時における応力−ひずみ関係
▲１▼ 応力σと初期降伏応力Ｙ₀ と等方硬化量Ｒと背応力αとの関係
前記式(1) において、ｘ軸を基準とする一軸応力状態を想定する。具体的には、複数の応力σのうち、σ_x はσとするとともに、それ以外のものはすべて０であると仮定し、偏差応力Ｓ_xy，Ｓ_yz，Ｓ_zxもすべて０であると仮定する。よって、偏差応力Ｓは、次の式(13)で表される。
【００６５】
【数２５】

【００６６】
さらに、前記式(4) より明らかなように、背応力増分ｄαは塑性ひずみ増分ｄε^P に比例し、また、ロイスの流動則に従い、塑性ひずみ増分ｄε^P の偏差成分は偏差応力Ｓに比例すると仮定する。よって、背応力αは、次の式(14)で表される。
【００６７】
【数２６】

【００６８】
それら式(13)および(14)で表された偏差応力Ｓおよび背応力αを前記式(1) に代入し、かつ、降伏条件として降伏関数ｆ＝０を与えると、次の式(15)が得られる。
【００６９】
【数２７】

【００７０】
▲２▼ 等方硬化量Ｒと相当塑性ひずみε^P _eqとの関係
前記式(10)を積分することにより、次の式(16)が得られる。
【００７１】
【数２８】

【００７２】
▲３▼ 係数Ｃ
前記係数Ｃは、引張時における相当塑性ひずみをε^P _eqで表すことにより、前記式(6) と同様に、次の式(17)で表される。
【００７３】
【数２９】

【００７４】
▲４▼ 第１背応力成分α₁
第１背応力成分α₁ は、前記式(3) に基づき、かつ、前記係数Ｃ’を用いることにより、次の式(18)で表される。
【００７５】
【数３０】

【００７６】
▲５▼ 第２背応力成分α₂
第２背応力成分α₂ は、相当塑性ひずみε^P _eqを変数とする関数により定義される。よって、第２背応力成分α₂ は、次の式(19)で定義される。
【００７７】
【数３１】

【００７８】
(4) 圧縮時における応力−ひずみ関係
▲１▼ 応力σと初期降伏応力Ｙ₀ と等方硬化量Ｒと背応力αとの関係
前記式(15)と同様にして、次の式(20)が得られる。
【００７９】
【数３２】

【００８０】
▲２▼ 等方硬化量Ｒと相当塑性ひずみε^P _eqとの関係
前記式(10)より、次の式(21)が得られる。
【００８１】
【数３３】

【００８２】
▲３▼ 係数Ｃ
係数Ｃは、圧縮時における相当塑性ひずみをε^P _eqc で表すことにより、前記式(6) と同様に、次の式(22)で表される。
【００８３】
【数３４】

【００８４】
▲４▼ 第１背応力成分α₁
第１背応力成分α₁ は、前記係数Ｃ’を用いることにより、次の式(23)で表される。
【００８５】
【数３５】

【００８６】
ただし、
α_1T：引張終了時における第１背応力成分α₁
【００８７】
▲５▼ 第２背応力成分α₂
第２背応力成分α₂ は、圧縮時における相当塑性ひずみε^P _eqc を変数とする関数により定義される。よって、第２背応力成分α₂ は、次の式(24)で定義される。
【００８８】
【数３６】

【００８９】
(5) 複合硬化モデルにおける各種材料定数および材料変数の決定
▲１▼ 等方硬化特性を定義する定数ｂと初期降伏応力Ｙ₀ の決定
前述のように、引張時には、次の式(25)が成立する。
【００９０】
【数３７】

【００９１】
また、引張時における降伏応力と圧縮時における降伏応力との差である応力差Δσを用いると、前記式(15)と(20)とから、次の式(26)が誘導される。
【００９２】
【数３８】

【００９３】
本実施形態においては、定数ｂと初期降伏応力Ｙ₀ の決定に先立ち、前述のように、試験片を引張方向に負荷して塑性変形させた後に除荷し、さらに、逆方向すなわち圧縮方向に負荷して塑性変形させるとともに、その間にその試験片の応力−ひずみ関係を実測することが行われる。すなわち、応力−ひずみ関係の実験値が存在するのである。よって、上記式(25)および(26)における定数および変数のうち、相当塑性ひずみε^P _eqと応力差Δσは既知であり、等方硬化量Ｒと定数ｂと初期降伏応力Ｙ₀ が未知である。一方、等方硬化量Ｒは、相当塑性ひずみε^P _eqと比例関係にあり、その比例定数がｂである。また、上記式(26)から明らかなように、等方硬化量Ｒと応力差Δσが判明すれば初期降伏応力Ｙ₀ を取得できる。
【００９４】
したがって、相当塑性ひずみε^P _eqに任意の値を離散的に与えるとともに、そのときの応力差Δσを実験値から取得し、一方、初期降伏応力Ｙ₀ は、相当塑性ひずみε^P _eqに依存しない固定値であるから、任意値に固定し、この状態で等方硬化量Ｒを上記式(26)により計算する。複数の相当塑性ひずみε^P _eqについて同様なことを繰返し、それにより、対応する複数の等方硬化量Ｒを計算する。そして、それら複数の相当塑性ひずみε^P _eqと複数の等方硬化量Ｒとに対して最小２乗法を適用することにより、比例定数ｂを計算する。比例定数ｂが計算されたならば、上記式(25)により、各相当塑性ひずみε^P _eqごとに、等方硬化量Ｒを計算し、さらに、上記式(26)により、初期降伏応力Ｙ₀ を計算する。
【００９５】
▲２▼ 移動硬化特性を定義する係数Ｈ，係数ａおよびＣ’の決定
本発明者らは、応力−ひずみ関係の実験値を検討することにより、第１背応力成分α₁ については、相当塑性ひずみε^P _eqが１％（０．０１）以下である極低ひずみ領域においては、相当塑性ひずみε^P _eqに応じて変化し、１％より大きい領域においては、相当塑性ひずみε^P _eqに応じて変化せず、一定となることに気がついた。また、第２背応力成分α₂ については、相当塑性ひずみε^P _eqが４％以下である低ひずみ領域においては、相当塑性ひずみε^P _eqに対する線形性が強く、４％より大きく、１５％以下である高い中ひずみ領域においては、非線形性が強く、１５％より大きい高ひずみ領域においては、相当塑性ひずみε^P _eqに対する線形性が強いことに気がついた。
【００９６】
要するに、図７に示すように、第１背応力成分α₁ については、相当塑性ひずみε^P _eqが境界値ε₀ 以下である極低ひずみ領域▲１▼と、境界値ε₀ より大きい非極低ひずみ領域▲２▼とに分けて検討することが適当であり、また、第２背応力成分α₂ については、相当塑性ひずみε^P _eqが境界値ε₁ 以下である低ひずみ領域▲３▼と、境界値ε₁ より大きく、境界値ε₂ 以下である中ひずみ領域▲４▼と、境界値ε₂ より大きい高ひずみ領域▲５▼とに分けて検討することが適当であることが判明したのである。
【００９７】
▲３▼ 引張時の移動硬化特性の決定
(i) 定数Ｈ₁ の決定
定数Ｈ₁ を決定することは、第２背応力成分α₂ を決定することを意味するが、第２背応力成分α₂ は、第１背応力成分α₁ と共同して背応力αを構成することから、定数Ｈ₁ を決定する際に、第１背応力成分α₁ が判明していることが必要である。
【００９８】
第１背応力成分α₁ は、極低ひずみ領域▲１▼においては、相当塑性ひずみε^P _eqに応じて変化するが、非極低ひずみ領域▲２▼においては一定値に収束する。すなわち、
α₁ ＝２ａ／３
となる。よって、定数Ｈ₁ を決定するに際し、非極低ひずみ領域▲２▼に属する低ひずみ領域▲３▼においては、前記式(2) より、
α₂ ＝α−２ａ／３
となる。また、背応力αは、前記式(15)から誘導されるように、
α＝２／３（σ−Ｙ₀ −Ｒ）
で表され、かつ、それら応力σと初期降伏応力Ｙ₀ と等方硬化量Ｒは前述のように、応力−ひずみ関係の実験値から、相当塑性ひずみε^P _eqの各値に関連付けて取得可能である。さらに、第２背応力成分α₂ は、前記式(4) に示すように、相当塑性ひずみε^P _eqと比例関係にある。それらを利用することにより、定数Ｈ₁ が決定される。具体的には、低ひずみ領域▲３▼に属する複数の離散値である相当塑性ひずみε^P _eqに対応する複数の背応力αが取得され、それら複数の相当塑性ひずみε^P _eqと複数の背応力αとに対して最小２乗法が適用されることにより、定数Ｈ₁ が決定される。
【００９９】
その結果、第２背応力成分α₂ が、次の式(27)で表される。
【０１００】
【数３９】

【０１０１】
(ii)係数ａとＣ’の決定
前記式(2) と(15)とから、次の式(28)が誘導される。
【０１０２】
【数４０】

【０１０３】
また、この式(28)における「α₁ 」を前記式(18)に代入することにより、次の式(29)が誘導される。
【０１０４】
【数４１】

【０１０５】
この式(29)は、次の式(30)に変形できる。
【０１０６】
【数４２】

【０１０７】
ここで、相当塑性ひずみε^P _eqが定数ε₀ より十分に大きいと仮定すれば、上記式(30)は、次の式(31)に変形できる。
【０１０８】
【数４３】

【０１０９】
ここに、「Ｃ’」は、基準ひずみε^P _T により決まり、相当塑性ひずみε^P _eqによって変化しないことから、定数として扱うことできる。また、「ａ」も、同様に、基準ひずみε^P _T により決まり、相当塑性ひずみε^P _eqによって変化しないことから、定数として扱うことできる。一方、上記式(31)は、相当塑性ひずみε^P _eqと、次の式(32)で表される項、すなわち、
【０１１０】
【数４４】

【０１１１】
とが比例関係にあることを示している。よって、相当塑性ひずみε^P _eqと、上記式(31)における垂直応力σ，初期降伏応力Ｙ₀ ，等方硬化量Ｒおよび第２背応力成分α₂ とに、前記複数の実験値のうち極低ひずみ領域▲１▼で取得されたものを代入し、相当塑性ひずみε^P _eqと、上記式(32)で表される項とが少ない誤差で比例関係が成立するように最小２乗法により、係数ａが決定される。このようにして係数ａが決定された後、前記式(29)を利用することにより、係数Ｃ’が決定される。
【０１１２】
このようにして係数ａおよび係数Ｃ’を決定することは、互いに異なる複数の基準ひずみε^P _T についてそれぞれ行われ、その結果、互いに対応する複数の係数ａと複数の基準ひずみε^P _T との関係と、互いに対応する複数の係数Ｃ’と複数の基準ひずみε^P _T との関係とが取得されることになる。
【０１１３】
一方、係数Ｃ’を決定するためには、前記式(11)を同定すること、すなわち、定数Ｃ’₁ およびＣ’₀ を決定することが必要である。そして、互いに対応する複数の係数Ｃ’と複数の基準ひずみε^P _T とに対して最小２乗法が適用されることにより定数Ｃ’₁ が前記式(11)における比例定数として決定され、さらに、その決定された定数Ｃ’₁ と、互いに対応する係数Ｃ’と基準ひずみε^P _T とを前記式(11)に代入することにより、定数Ｃ’₀ が決定される。
【０１１４】
また、同様に、係数ａを決定するためには、前記式(12)を同定すること、すなわち、定数ａ₁ およびａ₀ を決定することが必要である。そして、互いに対応する複数の係数ａと複数の基準ひずみε^P _T とに対して最小２乗法が適用されることにより定数ａ₁ が前記式(12)における比例定数として決定され、さらに、その決定された定数ａ₁ と、互いに対応する係数ａと基準ひずみε^P _T とを前記式(12)に代入することにより、定数ａ₀ が決定される。
【０１１５】
(iii) 定数ｐ，ｑの決定
前記式(4) と式(28)とから、次の式(33)が誘導される。
【０１１６】
【数４５】

【０１１７】
ただし、
ｙ：比例係数
【０１１８】
図８には、応力−ひずみ関係の実験値を用いて上記比例係数ｙがグラフで表されている。このグラフには複数の実験値に基づく比例係数ｙの複数の値がドットでプロットされるとともに、それらドットを近似する曲線が描かれている。前記式(8) が示すように、係数Ｈとｌｎ（ε^P _eq）との間に比例関係が成立するため、応力−ひずみ関係の実験値に基づき、最小２乗法を用いることにより、係数Ｈとｌｎ（ε^P _eq）との比例定数に相当する定数ｐが決定され、続いて、その決定された定数ｐと応力−ひずみ関係の実験値とを上記式(8) に代入することにより、定数ｑが決定される。
【０１１９】
前記式(8) のＨを前記式(4) に代入し、それにより得られた微分方程式を解くと、次の式(34)が得られる。
【０１２０】
【数４６】

【０１２１】
ここに、α₂ ¹は、前記式(27)に境界値ε₁ を代入して得られる定数である。
【０１２２】
(iv)定数Ｈ₂ の決定
定数Ｈ₂ は、図８のグラフのうち、相当塑性ひずみε^P _eqが境界値ε₂ 以上である高ひずみ領域▲５▼に対応する部分を参照することにより、決定される。
【０１２３】
前記式(9) のＨを前記式(4) に代入し、それにより得られた微分方程式を解くと、次の式(35)が得られる。
【０１２４】
【数４７】

【０１２５】
ここに、α₂ ²は、前記式(34)に境界値ε₂ を代入して得られる定数である。
【０１２６】
▲４▼ 圧縮時の移動硬化特性の決定
(i) 第１背応力成分α₁ の決定（係数ａおよびＣ’の決定）
前記式(2) と(20)とから、次の式(36)が誘導される。
【０１２７】
【数４８】

【０１２８】
この式(36)における「α₁ 」を前記式(23)に代入することにより、次の式(37)が誘導される。
【０１２９】
【数４９】

【０１３０】
この式(37)は、引張時に準じて、相当塑性ひずみε^P _eqと、次の式(38)で表される項、すなわち、
【０１３１】
【数５０】

【０１３２】
とが比例関係にあることを示している。よって、相当塑性ひずみε^P _eqと、上記式(36)における垂直応力σ，初期降伏応力Ｙ₀ ，等方硬化量Ｒおよび第２背応力成分α₂ とに、前記複数の実験値のうち極低ひずみ領域▲１▼で取得されたものを代入し、相当塑性ひずみε^P _eqと、上記式(38)で表される項とが少ない誤差で比例関係が成立するように最小２乗法により、係数ａが決定される。このようにして係数ａが決定された後、前記式(37)を利用することにより、係数Ｃ’が決定される。
【０１３３】
(ii)第２背応力成分α₂ の決定（係数Ｈの決定）
係数Ｈは、引張時と同じものを用いる。よって、第２背応力成分α₂ は、引張時における第２背応力成分α₂ との連続性を考慮し、低ひずみ領域▲３▼においては、次の式(39)で表される。
【０１３４】
【数５１】

【０１３５】
ただし、
α_2T：引張終了時における第２背応力成分α₂
【０１３６】
また、中ひずみ領域▲４▼においては、次の式(40)で表される。
【０１３７】
【数５２】

【０１３８】
また、高ひずみ領域▲５▼においては、次の式(41)で表される。
【０１３９】
【数５３】

【０１４０】
図９には、以上説明した同定を行うための前記塑性構成式同定プログラムがフローチャートで表されている。
【０１４１】
まず、ステップＳ１０１において、等方硬化特性を定義する定数ｂおよび初期降伏応力Ｙ₀ が決定される。次に、引張時における移動硬化特性を定義する係数Ｈ，係数ａおよび係数Ｃ’が決定される。その後、圧縮時における移動硬化特性を定義する係数Ｈ，係数ａおよび係数Ｃ’が決定される。以上で本プログラムの一回の実行が終了する。
【０１４２】
以上のようにして同定された前記式(1) の降伏関数ｆを用いることにより、材料の応力−ひずみ関係がシミュレートされる。
【０１４３】
図５および図６には、さらに、前記低ひずみ領域と高ひずみ領域とのそれぞれにつき、本実施形態による応力−ひずみ関係のシミュレート結果が実線グラフで示されている。それら２つのグラフから明らかなように、本実施形態によれば、そのシミュレートによる予測値が低ひずみ領域のみならず高ひずみ領域においても実験値に精度よく一致する。高ひずみ領域においてもバウシンガ効果が精度よくシミュレートされているのである。
【０１４４】
さらに、本実施形態においては、同定された降伏関数ｆと、成形シミュレーションにより予測された残存応力−ひずみ関係とに基づき、材料のプレス成形時におけるスプリングバック量が予測される。したがって、本実施形態によれば、高ひずみ領域においても精度よくバウシンガ効果を表現する降伏関数ｆを用いてスプリングバック量が予測されるため、その予測精度が向上する。
【０１４５】
以上の説明から明らかなように、本実施形態においては、図２のステップＳ１が「実験値取得工程」を構成し、ステップＳ２が「係数関数同定工程」および「関係同定工程」を構成し、ステップＳ３およびＳ４が互いに共同して「残存応力−ひずみ関係予測工程」を構成し、ステップＳ５およびＳ６が互いに共同して「スプリングバック量予測工程」を構成しているのである。また、前記式(1) が「降伏関数」であり、前記式(6) が、降伏関数ｆの２つの係数の一方である係数Ｃの「係数関数」であり、前記式(7) ないし(9) が、他方の係数Ｈの「係数関数」なのである。
【０１４６】
なお付言すれば、本実施形態によれば、応力−ひずみ関係を連続的にシミュレート可能な降伏関数ｆの基本的な形が予め定められていて、複数の離散的な実験値は主に、その降伏関数ｆにおける各種係数，定数を決定するために使用されるため、比較的少ない数の実験値から、応力−ひずみ関係を連続的にかつ高い精度でシミュレートすることが可能となる。
【０１４７】
さらに付言すれば、弾塑性体の挙動を有限要素法により数値解析するためのプログラムは一般に、テンソル量の増分で表現されている。したがって、材料硬化モデルをプレス成形等で用いられる有限要素法のプログラムに組み込むためには、材料硬化モデルもテンソル量の増分で表現されることが望ましい。これに対して、本実施形態においては、それの材料硬化モデルである複合硬化モデルがテンソル量の増分で表現されているため、材料硬化モデルを有限要素法のプログラムに容易に組み込むことが可能である。
【０１４８】
さらに付言すれば、本実施形態においては、「係数関数同定工程」が人間の介入を実質的に必要とせずにコンピュータにより自動的に実行されるようになっているが、人間が主体的に介入する一方、コンピュータを補助的に使用することによって実行することが可能である。
【０１４９】
以上、本発明の一実施形態を図面に基づいて詳細に説明したが、これは例示であり、前記〔発明が解決しようとする課題，課題解決手段および発明の効果〕の項に記載された態様を始めとして、当業者の知識に基づいて種々の変形，改良を施した形態で本発明を実施することが可能である。
【図面の簡単な説明】
【図１】本発明の一実施形態である応力−ひずみ関係シミュレート方法を含むスプリングバック量予測方法を実施するのに好適なスプリングバック量予測装置を示す系統図である。
【図２】上記スプリングバック量予測方法を示す工程図である。
【図３】図１のコンピュータ１２により実行される３つのプログラム間の関係を説明するためのブロック図である。
【図４】引張−圧縮時における応力σとひずみεとの関係を示すとともに、上記実施形態において用いられる複合硬化モデルにおける初期降伏応力Ｙ₀ と背応力αと等方硬化量Ｒとの関係を示すグラフである。
【図５】低ひずみ領域において材料の応力−ひずみ関係がシミュレートされる精度を、Chaboche-Rousselier により提案された一従来技術と上記実施形態とを実験値を参照しつつ対比するためのグラフである。
【図６】高ひずみ領域において材料の応力−ひずみ関係がシミュレートされる精度を、Chaboche-Rousselier により提案された一従来技術と上記実施形態とを実験値を参照しつつ対比するためのグラフである。
【図７】上記実施形態において第１背応力成分α₁ の特性が極低ひずみ領域▲１▼と非極低ひずみ領域▲２▼とに分けて解析されるとともに、第２背応力成分α₂ の特性が低ひずみ領域▲３▼と中ひずみ領域▲４▼と高ひずみ領域▲５▼とに分けて解析されることを示すグラフである。
【図８】上記実施形態において比例係数ｙが相当塑性ひずみε^P _eqに応じて変化する様子を示すグラフである。
【図９】図３における塑性構成式同定プログラムを示すフローチャートである。
【符号の説明】
１２ コンピュータ
１６ 外部記憶装置
３０ ＣＤ−ＲＯＭ
３２ ＦＤ[0001]
BACKGROUND OF THE INVENTION
  The present invention relates to a stress-strain relationship simulation method for simulating a stress-strain relationship of a material including the Bauschinger effect of the material., And a method for predicting the amount of springback using the simulation methodIt is about.
[0002]
[Prior art]
The following types of stress-strain relationship simulation methods already exist. It approximates the stress-strain relationship of a material with a composite hardening model that combines an isotropic hardening model and a kinematic hardening model, and uses the plastic strain of the material as a variable and has a coefficient of the plastic strain. By using the yield function that defines the yield surface of the material by using the back stress and the amount of isotropic hardening defined by the variables as variables, the stress-strain relationship of the material is continuously and This simulates the Bauschinger effect.
[0003]
In the yield function, the coefficient of plastic strain has conventionally been a constant. Furthermore, the relationship between the amount of isotropic hardening and the plastic strain is set so that the amount of isotropic hardening increases according to the plastic strain while the plastic strain is low, but does not increase when the plastic strain is high and converges to a constant value. It was.
[0004]
[Problems to be solved by the invention, problem-solving means, and effects of the invention]
For this reason, when this yield function is used, the material has a plastic strain greater than 0.02 (2%), for example, 0.15, 0.20. There was a problem that the effect could not be simulated accurately.
[0005]
  Against this background, the present invention has been made with the object of accurately simulating the Bauschinger effect of a material even in a high strain region, and the following aspects can be obtained by the present invention. As with the claims, each aspect is divided into sections, each section is numbered, and is described in a form that cites the numbers of other sections as necessary. This is to facilitate understanding of some of the technical features and combinations thereof described herein, and the technical features and combinations thereof described herein are limited to the following aspects. Should not be interpreted.
  Although some of the following items are no longer an embodiment of the present invention due to amendment of the claims, they are left as they are for the purpose of explaining the present invention.
[0006]
  (1) The stress-strain relationship of a material is approximated by a composite hardening model that combines an isotropic hardening model and a kinematic hardening model, and the back stress with the plastic strain of the material as a variable and the coefficient of the plastic strain. By using the yield function that defines the yield surface of the material by using the back stress defined by the function and the amount of isotropic hardening as variables, the stress-strain relationship of the material is continuously measured and the material A method for simulating the Bausinga effect of
  The material is loaded in one direction of the tensile direction and the compression direction and plastically deformed and then unloaded, and further, the material is plastically deformed by loading in the opposite direction, and during that time, the experimental value of the stress-strain relationship of the material Experimental value acquisition step of acquiring as a plurality of discrete values;
  The coefficient is a variable defined by a coefficient function having the plastic strain as a variable, and the coefficient function is identified based on a plurality of acquired experimental values.
  To simulate stress-strain relationshipLaw.
  According to this method, since the coefficient in the yield function is not a constant but a variable that is changed according to plastic strain, the Bauschinger effect can be accurately simulated even in a high strain region.
  Here, “plastic strain” may mean a plastic strain related to each coordinate axis in a multiaxial stress state, or a corresponding plastic strain when a multiaxial stress state is converted into a uniaxial stress state.
  The stress-strain relationship simulation method described in this section is based on the following: (a) The force that should be applied to the material during molding by changing the stress, strain, shape, etc. of the material during molding (including drawing). And the first step of simulating according to the stress-strain relationship of the material, and (b) the stress of the material based on the final values of the stress, strain, shape, etc. simulated by the first step The second step of predicting the amount of springback of the material according to the strain relationship is implemented in a mode applied to the second step of the springback amount prediction method, the mode applied in the first step, or the first step It is possible to implement in a mode that applies to both the second step and the second step.
  Further, the stress-strain relationship simulation method described in this section can be implemented so that the whole is independently executed by a human, or can be executed mainly by a computer. In particular, the “coefficient function identification step” may require human intervention or may be performed automatically by a computer without substantially requiring human intervention.
  The stress-strain relationship simulation method described in this section can be grasped as, for example, a stress-strain relationship analysis method or a stress-strain relationship definition method, or as a Bauschinger hardening analysis method.
  (2) The yield function f is
[Equation 8]

However,
  S_i(I = x, y, z): normal stress σ acting on a plane perpendicular to the coordinate axis corresponding to the index i_iDeviation stress on
  S_ij(I, j = x, y, z): Shear stress σ acting on the plane perpendicular to the coordinate axis corresponding to the index i in the direction of the coordinate axis corresponding to the index j_ijDeviation stress on
  α_i(I = x, y, z): Back stress acting on a plane perpendicular to the coordinate axis corresponding to the index i
  α_ij(I, j = x, y, z): Back stress acting on the plane perpendicular to the coordinate axis corresponding to the index i in the direction of the coordinate axis corresponding to the index j
  Y₀ : Initial yield stress, which is the yield stress when plastic deformation of material starts for the first time
  R: Isotropic amount of material
How to simulate the stress-strain relationship described in (1)Law.
  (3) A vector {α} representing the back stress α is
[Equation 9]

However,
  {Α₁ }: Vector representing the first back stress component which is a component depending on the nonlinearity of kinematic hardening in the back stress
  {Α₂ }: Vector representing the second back stress component which is a component depending on the linearity of kinematic hardening in the back stress
The first back stress component α₁ Increment dα₁ Representing the vector {dα₁ }But,
[Expression 10]

However,
  C: Coefficient representing the convergence speed of nonlinear kinematic hardening
  2a / 3: convergence value of nonlinear kinematic hardening
  {Dε^P}: Plastic strain increment dε^P _iA vector representing (i = x, y, z)
  ε^P _eq: Equivalent plastic strain
And the second back stress component α₂ Increment dα₂ Representing the vector {dα₂ }But,
## EQU11 ##

However,
  H: Coefficient representing the magnitude of linear kinematic hardening
How to simulate the stress-strain relationship described in (2)Law.
  (4) Further, the relationship between the amount of isotropic hardening and the plastic strain is such that the amount of isotropic hardening increases according to the plastic strain regardless of the magnitude of the plastic strain, and the amount of isotropic hardening and the plastic strain. The stress-strain relationship simulation method according to any one of (1) to (3), including a relationship identification step for identifying a relationship with strain based on a plurality of acquired experimental values.[Part of this item corresponds to claim 3].
  According to this method, the amount of isotropic hardening increases according to the plastic strain regardless of the magnitude of the plastic strain. Therefore, the amount of isotropic hardening changes depending on the plastic strain even in the high strain region, and the Bausinger effect is accurately obtained. Convenient to simulate.
  (5) The increment dR of the isotropic hardening amount R is
[Expression 12]

However,
  b: Constant indicating the magnitude of isotropic hardening
The stress-strain relationship simulation method described in item (4)[Part of this item corresponds to claim 4].
  (6) The stress-strain relationship simulation method according to any one of (3) to (5), wherein the coefficient C is defined by a function having a tendency to decrease according to the plastic strain.[In this section (3) From the subordinate to the term (2) Terms and (3) Excluding the limitation of the specific formula of the term corresponds to a part of claim 1).
  (7) The function representing the coefficient C is equivalent plastic strain ε as the plastic strain.^P _eqIs equivalent plastic strain ε^P _eq1 / √ (ε^P _eq(6) The stress-strain relationship simulation method according to (6).
  (8) The coefficient C is
[Formula 13]

However,
  ε^P _eq: Equivalent plastic strain
  C ', ε₀ :constant
The stress-strain relationship simulation method according to any one of (3) to (7)[Corresponds to a part of claim 2].
  (9) The stress-strain relationship simulation method according to any one of (3) to (8), wherein the coefficient H is defined by a function having a tendency to increase according to the plastic strain.[In this section (3) From the subordinate to the term (2) Terms and (3) Excluding the limitation of the specific formula of the term corresponds to a part of claim 1).
  (10) The function representing the coefficient H is the equivalent plastic strain ε as the plastic strain.^P _eqIs equivalent plastic strain ε^P _eqIs the natural logarithm ln (ε^P _eq) Or common logarithm log (ε^P _eqThe stress-strain relationship simulation method according to item (9), which is incorporated in the form of
  (11) The coefficient H is
[Expression 14]

However,
  ε₁ :Boundary value
  ε₂ : Boundary value (> ε₁ )
  H₁ , P, q, H₂ :constant
The stress-strain relationship simulation method according to any one of (3) to (10)[Corresponds to a part of claim 2].
  (12) A recording medium on which a computer-readable program is recorded so as to execute the stress-strain relationship simulation method according to any one of (1) to (11).
  Here, “recording medium” includes, for example,flexibleThere are disks, magnetic tapes, magnetic disks, magnetic drums, magnetic cards, optical disks, magneto-optical disks, ROMs, CD-ROMs, IC cards, punched tapes, and the like. This description also applies to a “recording medium” described later.
  (13) By using the yield function according to any one of the items (1) to (11), a method for predicting the amount of springback, which is the amount by which the material elastically recovers after the end of plastic working,
  A specimen having substantially the same material as the material is loaded in one direction in the tensile direction and the compression direction and plastically deformed, then unloaded, and further loaded in the opposite direction to be plastically deformed. An experimental value acquisition step of acquiring the experimental value of the stress-strain relationship of the test piece as a plurality of discrete values;
  The coefficient of the back stress in the yield function is a variable defined by a coefficient function having the plastic strain as a variable, and the coefficient function is identified based on a plurality of acquired experimental values. Process,
  A residual stress-strain relationship prediction step for predicting a stress-strain relationship remaining in the material at the end of the plastic working;
  A springback amount prediction step of predicting the springback amount based on the predicted residual stress-strain relationship and the yield function including the coefficient defined by the identified coefficient function;
  A method for predicting the amount of springback characterized in that[Part of this item corresponds to claim 5].
  In plastic working such as press forming, the material is unloaded after being loaded, but after unloading, the material has springback, and the shape at the end of unloading of the material (final shape) It will be different from the shape at the time of loading. Therefore, in order to ensure the accuracy of the shape of the material after plastic working, the shape of the tool for plastic working, such as a press die, is modified to a tool shape that allows for the amount of springback in the tool shape corresponding to the target shape of the material. It is necessary to do. This tool correction can be realized by trying plastic processing using a tool whose shape is provisional and actually measuring the springback amount of the material, and correcting the tool shape in consideration of the actual measurement value. However, if the amount of springback can be predicted with high accuracy without trying plastic working, the time required for tool correction is shortened. If the time required for tool correction is shortened, the production preparation lead time is shortened and the production efficiency is improved. However, conventionally, it has been difficult to accurately predict the amount of springback of the material. As described above, this is because the conventional method cannot accurately simulate the Bauschinger effect of a material in a high strain region.
  On the other hand, according to the springback amount prediction method described in this section, it is possible to accurately simulate the bauschinger effect of the material even in a high strain region. On the other hand, if the simulation accuracy of the Bauschinger effect of the material is increased, the prediction accuracy of the springback amount is also increased. Therefore, according to this method, the springback amount can be accurately predicted.
  (14) Further, the relationship between the amount of isotropic hardening and the plastic strain is such that the amount of isotropic hardening increases according to the plastic strain regardless of the magnitude of the plastic strain, and the amount of isotropic hardening and the plastic strain. Including a relationship identification step for identifying a relationship with strain based on a plurality of acquired experimental values, and the step of predicting the amount of springback includes: (a) a predicted residual stress-strain relationship; and (b) The stress according to (13), wherein the springback amount is predicted based on the yield function including the coefficient defined by the identified coefficient function and a coefficient representing the identified relationship. Strain relationship simulation method.
  According to this method, the amount of isotropic hardening increases according to the plastic strain regardless of the magnitude of the plastic strain. Therefore, the amount of isotropic hardening changes depending on the plastic strain even in the high strain region, and the Bausinger effect is accurately obtained. Convenient to simulate.
  (15) A recording medium on which a computer-executable program for recording the springback amount predicting method according to (13) or (14) is recorded in a computer-readable manner[Part of this item corresponds to claim 6].
[0007]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, one of more specific embodiments of the present invention will be described in detail with reference to the drawings.
[0008]
One embodiment of the present invention is a springback amount prediction method including a stress-strain relationship simulation method. FIG. 1 shows a springback amount prediction apparatus (hereinafter referred to as a springback amount prediction device) suitable for implementing the springback amount prediction method. Simply referred to as “prediction device”).
[0009]
The prediction device includes an input device 10, a computer 12, an output device 14, and an external storage device 16. The input device 10 is configured to include a mouse, a keyboard, and the like. The computer 12 includes a processor 20 such as a CPU, a memory 22 such as a ROM, a RAM, and a hard disk, and a bus 24 that connects the processor 20 and the memory 22. The output device 14 is configured to include a display, a printer, a plotter, and the like. The external storage device 16 can be loaded with a recording medium such as a CD-ROM 30 and a writable FD 32 (flexible disk). In the loaded state, reading and writing of data with respect to the recording medium are performed as necessary. .
[0010]
In the present embodiment, a program necessary for predicting the springback amount is stored in the CD-ROM 30, and data necessary for the prediction is stored in the FD 32. Programs necessary for predicting the springback amount include a plastic constitutive equation identification program, a forming simulation program, and a springback amount prediction program. When the springback amount is predicted, necessary programs and data are read from the CD-ROM 30 and the FD 32, transferred to the RAM or hard disk of the computer 12, and then the program is executed by the processor 20.
[0011]
FIG. 2 is a process diagram illustrating the springback amount prediction method according to the present embodiment.
[0012]
If it demonstrates roughly, the springback amount prediction method will predict the amount of springback which arises in the material, after pressing the material which is a steel plate. The press molding is performed so that a portion having a relatively high plastic strain is generated, such as a plastic strain of 0.15 and 0.20. In addition, the stress-strain relationship simulation method included in the springback amount prediction method uses a test piece that is substantially the same material as the material to be press-molded, so that the stress-strain relationship test of the test piece is performed. A value is acquired as a plurality of discrete values, and a plastic constitutive equation, which is an equation constituting the stress-strain relationship of the test piece, is identified based on the acquired plurality of experimental values. This plastic constitutive equation describes the stress-strain relationship continuously. In addition, the springback amount prediction method uses a finite element method to perform a molding simulation that simulates the press molding of the material in order to predict the stress-strain relationship remaining in the material immediately after the press molding, The springback amount of the material is predicted by using the finite element method and the identified plastic constitutive equation.
[0013]
More specifically, in the present embodiment, as shown in the figure, first, in step S1, a test specimen having substantially the same material as the material to be press-molded is prepared by the operator, Thereafter, by using a testing machine, a tensile-compression test is performed on the test piece in a uniaxial stress state. Specifically, this test involves unloading the specimen after it is plastically deformed by loading it in the tensile direction, and then plastically deforming by loading it in the opposite direction, that is, in the compression direction. Experimental values of the stress-strain relationship are obtained as a plurality of discrete values.
[0014]
Next, in step S2, the plastic constitutive equation identification program is executed by the computer 12, thereby identifying the plastic constitutive equation based on the plurality of acquired experimental values. This identification will be described in detail later.
[0015]
Thereafter, in step S3, data necessary for executing the molding simulation, that is, data for numerical analysis is created by the operator. As shown in FIG. 3, the numerical analysis data includes mesh data for dividing the material into a plurality of meshes, conditions relating to stress applied to the material during press molding, and boundary conditions relating to the plurality of meshes. It is out.
[0016]
Subsequently, in step S4, the computer 12 executes the forming simulation program, thereby predicting the stress-strain relationship remaining in the material immediately after press forming. Furthermore, the shape of the material immediately after press molding is also predicted. An example of this molding simulation program is manufactured by LSTC and sold by Japan Research Institute under the name “LS-DYNA3D”.
[0017]
Thereafter, in step S5, the stress-strain relationship of the material whose springback amount is to be predicted is defined in the springback amount prediction program based on the plastic constitutive equation identified in step S2. This springback amount prediction program is designed so that the stress-strain relationship of the material can be arbitrarily defined.
[0018]
Subsequently, in step S6, the computer 12 executes the above-described springback amount prediction program, thereby predicting the springback amount. An example of this springback amount prediction program is manufactured by the Japan Research Institute and sold under the name “JOH-NIKE3D” by the Institute.
[0019]
This completes one execution of the springback amount prediction method. FIG. 3 conceptually shows the relationship between these three programs.
[0020]
Here, the identification of the plastic constitutive equation will be described in detail.
[0021]
1. Conventional composite curing model
The journal "Journal of Pressure Vessel Technology (MAY 1983, Vol. 105 / 154-159)" published a paper entitled "On the Plastic and Viscoplastic Constitutive Equations" by JLChaboche and G. Rousselier. ing. In this paper, a plastic constitutive equation, which is a constitutive equation that gives the stress-strain relationship of an elastoplastic body, is described. This plastic constitutive equation can express the Bauschinger effect of a material by approximating the stress-strain relationship of the material with a composite curing model in which an isotropic hardening model and a kinematic hardening model are combined.
[0022]
Here, “isotropic hardening” means that the yield surface is expanded isotropically without changing its center position and shape due to work hardening. “Moving hardening” means that the yield curved surface moves without changing its shape and size due to work hardening. The “yield surface” refers to a curved surface (including a polygon) representing the yield function in the stress space or stress plane, and the “yield function” refers to the stress component σ when the material yields._ijA function of (i, j = x, y, z).
[0023]
The yield function f as a plastic constitutive equation proposed by Chaboche-Rousselier is expressed by the following equation (1).
[0024]
[Expression 15]

[0025]
However,
S_i(I = x, y, z): normal stress σ acting on a plane perpendicular to the coordinate axis corresponding to the index i_iDeviation stress on
S_ij(I, j = x, y, z): Shear stress S acting on the plane perpendicular to the coordinate axis corresponding to the index i in the direction of the coordinate axis corresponding to the index j_ij(= Τ_ij) Deviation stress α_i(I = x, y, z): Back stress acting on a plane perpendicular to the coordinate axis corresponding to the index i
α_ij(I, j = x, y, z): Back stress acting on the plane perpendicular to the coordinate axis corresponding to the index i in the direction of the coordinate axis corresponding to the index j
Y₀: Initial yield stress, which is the yield stress when plastic deformation of material starts for the first time
R: Isotropic amount of material
[0026]
The above formula (1) is established on the premise of the Mises-type yield condition, isotropic hardening and linear hardening.
[0027]
In this equation (1), the back stress α_i, Α_ijRepresents the amount of movement of the center of the yield surface, and can be said to be the amount of kinematic hardening in the kinematic hardening model. In FIG. 4, the stress-strain relationship during tension-compression is represented by a graph, and the back stress α, the isotropic hardening amount R, and the initial yield stress Y₀The relationship is shown. That is, in this figure, the composite curing model is conceptually shown as a graph.
[0028]
In the uniaxial stress state, the back stress α is obtained by Chaboche-Rousselier._iA vector {α} representing (i = x, y, z) is defined by the following equation (2).
[0029]
[Expression 16]

[0030]
However,
{Α₁}: Vector representing the first back stress component which is a component depending on the nonlinearity of kinematic hardening in the back stress
{Α₂}: Vector representing the second back stress component which is a component depending on the linearity of kinematic hardening in the back stress
[0031]
Here, the first back stress component α₁Vector {α₁} Is defined by the following equation (3).
[0032]
[Expression 17]

[0033]
However,
C: Coefficient (constant) representing the speed of convergence of nonlinear kinematic hardening
2a / 3: convergence value of non-linear kinematic hardening (constant)
{Dε^P}: Plastic strain increment dε^P _iA vector representing (i = x, y, z)
dε^P _eq: Equivalent plastic strain increment
[0034]
The second back stress component α₂Vector {α₂} Is defined by the following equation (4).
[0035]
[Expression 18]

[0036]
However,
H: Coefficient (constant) indicating the magnitude of linear kinematic hardening
[0037]
Further, the increment dR of the isotropic hardening amount R is defined by the following equation (5).
[0038]
[Equation 19]

[0039]
However,
b: Coefficient (constant) indicating the magnitude of isotropic hardening
Q: Convergence value of isotropic hardening (constant)
[0040]
2. Problems of the conventional composite curing model
FIG. 5 shows the plastic strain ε^PThe stress-strain relationship is shown in the case where the material is pulled and plastically deformed and then unloaded, and the material is further compressed and plastically deformed in order to shift to a low strain region of less than 0.02. The experimental value of the stress-strain relationship is represented by a graph in which a plurality of white circles are arranged, while the predicted value based on the conventional composite hardening model proposed by Chaboche-Rousselier is represented by a broken line graph. As is clear from these graphs, in the low strain region, the predicted value based on the conventional composite hardening model matches the experimental value with high accuracy, and the Bauschinger effect can be simulated with high accuracy.
[0041]
In FIG. 6, in order to shift to a high strain region where the plastic strain is larger than 0.02 and smaller than 0.154, the material is pulled and plastically deformed, then unloaded, and the material is further compressed and plastically deformed. The stress-strain relationship is shown. The experimental value is represented by a graph in which a plurality of white circles are arranged, while the predicted value based on the conventional composite curing model by Chaboche-Rousselier is represented by a broken line graph. As is apparent from these graphs, in the high strain region, the predicted value based on the conventional composite hardening model does not coincide with the experimental value with high accuracy, and the Bauschinger effect cannot be accurately simulated.
[0042]
3. Composite curing model in this embodiment
(1) Model design goals
The Bauschinger effect can be accurately simulated not only in the low strain region but also in the high strain region.
[0043]
(2) Overview of composite curing model
(1) Variable of coefficient C
The coefficient C in the equation (3) is not a constant but equivalent plastic strain ε^P _eqIs a variable defined by a function (which is an example of a coefficient function). Coefficient C is equivalent plastic strain ε^P _eqThere is a tendency to decrease according to. Specifically, the coefficient C is expressed by the following equation (6).
[0044]
[Expression 20]

[0045]
However,
C ′: coefficient identified using experimental values of stress-strain relationship (described later)
ε₀: Equivalent plastic strain ε^P _eqIs a constant for preventing the denominator of the above formula (6) from becoming 0 when 10 is 0 (for example, 10^-6Degree value)
[0046]
(2) Variable of coefficient H
The coefficient H in the equation (4) is not a constant but an equivalent plastic strain ε^P _eqIs a variable defined by a function (which is an example of a coefficient function). Coefficient H is equivalent plastic strain ε^P _eqThere is a tendency to increase in accordance with. Specifically, the coefficient H is expressed by the following equations (7) to (9).
[0047]
[Expression 21]

[0048]
However,
ε₁:Boundary value
H₁p, q, H₂: Constants identified using experimental values of stress-strain relationship
ε₂: Boundary value (> ε₁)
[0049]
(3) Deformation of constitutive equation for isotropic hardening amount R
Equivalent plastic strain ε^P _eqIsotropic hardening amount R is equivalent plastic strain ε regardless of the size^P _eqAccordingly, the amount of isotropic hardening R is equivalent to the plastic strain ε not only in the low strain region but also in the high strain region.^P _eqTo increase according to. This is convenient for accurately simulating the Bauschinger effect even in a high strain region. Specifically, the equation (5) is transformed into the following equation (10).
[0050]
[Expression 22]

[0051]
As a result, the isotropic hardening amount R and the equivalent plastic strain ε^P _eqThe equivalent plastic strain ε^P _eqA proportional relationship is established regardless of the size of.
[0052]
(5) Determination of coefficient C '
The coefficient C ′ in the equation (6) is not a constant, but equivalent plastic strain ε at the end of tension (at the start of unloading).^P _eqIs the reference strain ε^P _TIs a variable defined by a function with. Therefore, the coefficient C ′ is uniquely determined for each tension-compression state. In the same tension-compression state, the equivalent plastic strain ε^P _eqDoes not change even if changes.
[0053]
A function representing the coefficient C ′ is represented by the following equation (11), for example.
[0054]
[Expression 23]

[0055]
However,
C ’₁, C ’₀: Constants identified using experimental values of stress-strain relationship
[0056]
The coefficient C ′ is the reference strain ε during tension.^P _TIs 0, the coefficient C ′₀Is equal to
[0057]
Further, the coefficient C ′ can be defined by another formula, for example, an n (≧ 2) order reference strain ε.^P _TIt can be defined by a polynomial including the terms.
[0058]
(6) Determination of coefficient a
The coefficient a in the equation (3) is not a constant, but a reference strain ε^P _TIs a variable defined by a function with. Therefore, the coefficient a is uniquely determined for each tension-compression state, similarly to the coefficient C ′, and in the same tension-compression state, the equivalent plastic strain ε^P _eqDoes not change even if changes.
[0059]
The function representing the coefficient a is represented by the following equation (12), for example.
[0060]
[Expression 24]

[0061]
However,
a₁, A₀: Constants identified using experimental values of stress-strain relationship
[0062]
The coefficient a is the reference strain ε during tension.^P _TIs 0, the coefficient a₀Is equal to
[0063]
Further, the coefficient a can be defined by another formula, for example, n (≧ 2) th order reference strain ε^P _TIt can be defined by a polynomial including the terms.
[0064]
(3) Stress-strain relationship during tension
(1) Stress σ and initial yield stress Y₀Between isotropic hardening amount R and back stress α
In the equation (1), a uniaxial stress state with the x axis as a reference is assumed. Specifically, among a plurality of stresses σ,_xIs assumed to be σ and all others are zero, and the deviation stress S_xy, S_yz, S_zxAre all zero. Therefore, the deviation stress S is expressed by the following equation (13).
[0065]
[Expression 25]

[0066]
Further, as apparent from the above equation (4), the back stress increment dα is the plastic strain increment dε.^PAnd in accordance with Royce's flow law, the plastic strain increment dε^PIt is assumed that the deviation component of is proportional to the deviation stress S. Therefore, the back stress α is expressed by the following equation (14).
[0067]
[Equation 26]

[0068]
When the deviation stress S and the back stress α expressed by the equations (13) and (14) are substituted into the equation (1) and the yield function f = 0 is given as the yield condition, the following equation (15) Is obtained.
[0069]
[Expression 27]

[0070]
(2) Isotropic hardening amount R and equivalent plastic strain ε^P _eqRelationship with
By integrating the equation (10), the following equation (16) is obtained.
[0071]
[Expression 28]

[0072]
(3) Coefficient C
The coefficient C represents the equivalent plastic strain during tension ε^P _eqIs represented by the following formula (17), similarly to the formula (6).
[0073]
[Expression 29]

[0074]
(4) First back stress component α₁
First back stress component α₁Is expressed by the following equation (18) based on the equation (3) and using the coefficient C ′.
[0075]
[30]

[0076]
(5) Second back stress component α₂
Second back stress component α₂Is equivalent plastic strain ε^P _eqIs defined by a function with Therefore, the second back stress component α₂Is defined by the following equation (19).
[0077]
[31]

[0078]
(4) Stress-strain relationship during compression
(1) Stress σ and initial yield stress Y₀Between isotropic hardening amount R and back stress α
The following equation (20) is obtained in the same manner as the equation (15).
[0079]
[Expression 32]

[0080]
(2) Isotropic hardening amount R and equivalent plastic strain ε^P _eqRelationship with
From the equation (10), the following equation (21) is obtained.
[0081]
[Expression 33]

[0082]
(3) Coefficient C
Coefficient C is the equivalent plastic strain during compression ε^P _eqcIn the same manner as the above formula (6), it is represented by the following formula (22).
[0083]
[Expression 34]

[0084]
(4) First back stress component α₁
First back stress component α₁Is expressed by the following equation (23) by using the coefficient C ′.
[0085]
[Expression 35]

[0086]
However,
α_1T: First back stress component α at the end of tension₁
[0087]
(5) Second back stress component α₂
Second back stress component α₂Is the equivalent plastic strain ε during compression^P _eqcIs defined by a function with Therefore, the second back stress component α₂Is defined by the following equation (24).
[0088]
[Expression 36]

[0089]
(5) Determination of various material constants and material variables in the composite curing model
(1) Constant b that defines isotropic hardening characteristics and initial yield stress Y₀Decision
As described above, the following equation (25) is established at the time of tension.
[0090]
[Expression 37]

[0091]
Further, when the stress difference Δσ, which is the difference between the yield stress during tension and the yield stress during compression, is used, the following equation (26) is derived from the equations (15) and (20).
[0092]
[Formula 38]

[0093]
In the present embodiment, the constant b and the initial yield stress Y₀Prior to the determination of the test piece, as described above, the specimen is loaded in the tensile direction and plastically deformed and then unloaded, and further in the opposite direction, that is, in the compression direction, plastically deformed. The stress-strain relationship is actually measured. That is, there is an experimental value of the stress-strain relationship. Therefore, among the constants and variables in the above equations (25) and (26), the equivalent plastic strain ε^P _eqAnd stress difference Δσ are known, isotropic hardening amount R, constant b, and initial yield stress Y₀Is unknown. On the other hand, isotropic hardening amount R is equivalent plastic strain ε^P _eqThe proportionality constant is b. Further, as apparent from the above equation (26), if the isotropic hardening amount R and the stress difference Δσ are found, the initial yield stress Y₀Can be obtained.
[0094]
Therefore, equivalent plastic strain ε^P _eqIs given discretely, and the stress difference Δσ at that time is obtained from experimental values, while the initial yield stress Y₀Is equivalent plastic strain ε^P _eqTherefore, the isotropic hardening amount R is calculated by the above equation (26). Multiple equivalent plastic strains ε^P _eqThe same is repeated for, thereby calculating a corresponding plurality of isotropic cure amounts R. And these multiple equivalent plastic strains ε^P _eqThe proportionality constant b is calculated by applying the least square method to the isotropic hardening amount R. Once the proportionality constant b is calculated, each equivalent plastic strain ε can be calculated by the above equation (25).^P _eqFor each, the isotropic hardening amount R is calculated, and the initial yield stress Y is calculated by the above equation (26).₀Calculate
[0095]
(2) Determination of coefficient H, coefficient a, and C 'that define kinematic hardening characteristics
By examining experimental values of the stress-strain relationship, the present inventors have obtained the first back stress component α.₁Is equivalent plastic strain ε^P _eqIn the extremely low strain region where the ratio is 1% (0.01) or less, the equivalent plastic strain ε^P _eqIn a region greater than 1%, the equivalent plastic strain ε^P _eqI noticed that it did not change according to the situation, and was constant. The second back stress component α₂Is equivalent plastic strain ε^P _eqIn the low strain region where the ratio is 4% or less, the equivalent plastic strain ε^P _eqIn the high medium strain region, which is greater than 4% and 15% or less, the non-linearity is strong, and in the high strain region greater than 15%, the equivalent plastic strain ε^P _eqI noticed that the linearity is strong.
[0096]
In short, as shown in FIG. 7, the first back stress component α₁Is equivalent plastic strain ε^P _eqIs the boundary value ε₀Extremely low strain region (1) below and boundary value ε₀It is appropriate to consider separately the larger non-very low strain region (2), and the second back stress component α₂Is equivalent plastic strain ε^P _eqIs the boundary value ε₁The following low strain region (3) and boundary value ε₁Larger, boundary value ε₂The medium strain region (4) below and the boundary value ε₂It has been found that it is appropriate to study separately in a larger high strain region (5).
[0097]
(3) Determination of kinematic hardening characteristics during tension
(i) Constant H₁Decision
Constant H₁Determining the second back stress component α₂The second back stress component α₂Is the first back stress component α₁Since the back stress α is composed in conjunction with₁In determining the first back stress component α₁Need to be known.
[0098]
First back stress component α₁Is equivalent plastic strain ε in the extremely low strain region (1).^P _eqHowever, it converges to a constant value in the non-very low strain region (2). That is,
α₁= 2a / 3
It becomes. Therefore, the constant H₁In the low strain region {circle around (3)} belonging to the non-very low strain region {circle around (2)},
α₂= Α-2a / 3
It becomes. Further, the back stress α is derived from the equation (15),
α = 2/3 (σ−Y₀-R)
And the stress σ and the initial yield stress Y₀The isotropic hardening amount R is equivalent to the equivalent plastic strain ε from the experimental value of the stress-strain relationship as described above.^P _eqCan be obtained in association with each value of. Further, the second back stress component α₂Is equivalent plastic strain ε as shown in the equation (4).^P _eqIs proportional to By using them, the constant H₁Is determined. Specifically, the equivalent plastic strain ε, which is a plurality of discrete values belonging to the low strain region (3).^P _eqA plurality of back stress α corresponding to the^P _eqAnd a plurality of back stresses α, the least square method is applied, so that the constant H₁Is determined.
[0099]
As a result, the second back stress component α₂Is expressed by the following equation (27).
[0100]
[39]

[0101]
(ii) Determination of coefficients a and C '
From the equations (2) and (15), the following equation (28) is derived.
[0102]
[Formula 40]

[0103]
Further, in this equation (28), “α₁Is substituted into the equation (18), the following equation (29) is derived.
[0104]
[Expression 41]

[0105]
This equation (29) can be transformed into the following equation (30).
[0106]
[Expression 42]

[0107]
Where the equivalent plastic strain ε^P _eqIs a constant ε₀Assuming that it is sufficiently larger, the above equation (30) can be transformed into the following equation (31).
[0108]
[Expression 43]

[0109]
Where “C ′” is the reference strain ε^P _TEquivalent plastic strain ε^P _eqCan be treated as a constant. Similarly, “a” is the reference strain ε.^P _TEquivalent plastic strain ε^P _eqCan be treated as a constant. On the other hand, the above equation (31) is equivalent plastic strain ε^P _eqAnd a term represented by the following equation (32), that is,
[0110]
(44)

[0111]
And are in a proportional relationship. Therefore, equivalent plastic strain ε^P _eqAnd the normal stress σ and the initial yield stress Y in the above equation (31)₀, Isotropic hardening amount R and second back stress component α₂And substituting the value obtained in the extremely low strain region {circle around (1)} among the plurality of experimental values.^P _eqThe coefficient a is determined by the method of least squares so that the proportional relationship is established with a small error between the term represented by the above equation (32). After the coefficient a is determined in this way, the coefficient C ′ is determined by using the equation (29).
[0112]
Determining the coefficient a and the coefficient C ′ in this way means that a plurality of different reference strains ε^P _TAs a result, a plurality of coefficients a and a plurality of reference strains ε corresponding to each other are obtained.^P _TAnd a plurality of coefficients C ′ and a plurality of reference strains ε corresponding to each other.^P _TAnd the relationship is acquired.
[0113]
On the other hand, in order to determine the coefficient C ′, the above equation (11) is identified, that is, the constant C ′.₁And C '₀It is necessary to determine A plurality of coefficients C ′ and a plurality of reference strains ε corresponding to each other.^P _TAnd the constant C ′₁Is determined as the proportionality constant in the equation (11), and the determined constant C ′₁And the corresponding coefficient C ′ and the reference strain ε^P _TIs substituted into the equation (11) to obtain a constant C ′₀Is determined.
[0114]
Similarly, in order to determine the coefficient a, the above equation (12) is identified, that is, the constant a₁And a₀It is necessary to determine A plurality of coefficients a and a plurality of reference strains ε corresponding to each other^P _TIs applied to the constant a₁Is determined as the proportionality constant in the equation (12), and the determined constant a₁And the corresponding coefficient a and the reference strain ε^P _TIs substituted into the equation (12) to obtain a constant a₀Is determined.
[0115]
(iii) Determination of constants p and q
From the equations (4) and (28), the following equation (33) is derived.
[0116]
[Equation 45]

[0117]
However,
y: Proportional coefficient
[0118]
In FIG. 8, the proportionality coefficient y is represented by a graph using experimental values of the stress-strain relationship. In this graph, a plurality of values of the proportionality coefficient y based on a plurality of experimental values are plotted with dots, and a curve approximating these dots is drawn. As the above equation (8) shows, the coefficients H and ln (ε^P _eq)), The coefficient H and ln (ε) are obtained by using the least square method based on the experimental value of the stress-strain relationship.^P _eqThe constant p corresponding to the proportionality constant is determined, and then the constant q is determined by substituting the determined constant p and the experimental value of the stress-strain relationship into the above equation (8). .
[0119]
Substituting H in the equation (8) into the equation (4) and solving the resulting differential equation, the following equation (34) is obtained.
[0120]
[Equation 46]

[0121]
Where α₂ ¹Is the boundary value ε in equation (27).₁Is a constant obtained by substituting
[0122]
(iv) Constant H₂Decision
Constant H₂Is equivalent plastic strain ε in the graph of FIG.^P _eqIs the boundary value ε₂This is determined by referring to the portion corresponding to the high strain region (5).
[0123]
Substituting H in the equation (9) into the equation (4) and solving the differential equation thus obtained, the following equation (35) is obtained.
[0124]
[Equation 47]

[0125]
Where α₂ ²Is the boundary value ε in the equation (34).₂Is a constant obtained by substituting
[0126]
(4) Determination of kinematic hardening characteristics during compression
(i) First back stress component α₁(Determination of coefficients a and C ')
From the equations (2) and (20), the following equation (36) is derived.
[0127]
[Formula 48]

[0128]
In this equation (36),₁Is substituted into the equation (23), the following equation (37) is derived.
[0129]
[Equation 49]

[0130]
This equation (37) is equivalent to the equivalent plastic strain ε^P _eqAnd a term represented by the following equation (38), that is,
[0131]
[Equation 50]

[0132]
And are in a proportional relationship. Therefore, equivalent plastic strain ε^P _eqAnd the normal stress σ and the initial yield stress Y in the above equation (36)₀, Isotropic hardening amount R and second back stress component α₂And substituting the value obtained in the extremely low strain region {circle around (1)} among the plurality of experimental values.^P _eqThe coefficient a is determined by the least square method so that the proportional relationship is established with a small error between the term represented by the above equation (38). After the coefficient a is determined in this way, the coefficient C ′ is determined by using the equation (37).
[0133]
(ii) Second back stress component α₂Determination (determination of coefficient H)
The coefficient H is the same as that used during tension. Therefore, the second back stress component α₂Is the second back stress component α during tension₂In the low strain region (3), the following equation (39) is used.
[0134]
[Equation 51]

[0135]
However,
α_2T: Second back stress component α at the end of tension₂
[0136]
Further, the medium strain region (4) is expressed by the following equation (40).
[0137]
[Formula 52]

[0138]
In the high strain region (5), it is expressed by the following equation (41).
[0139]
[Equation 53]

[0140]
FIG. 9 is a flowchart showing the plastic constitutive equation identification program for performing the identification described above.
[0141]
First, in step S101, a constant b and initial yield stress Y that define isotropic hardening characteristics₀Is determined. Next, a coefficient H, a coefficient a, and a coefficient C ′ that define kinematic hardening characteristics during tension are determined. Thereafter, a coefficient H, a coefficient a, and a coefficient C ′ that define kinematic hardening characteristics during compression are determined. This completes one execution of this program.
[0142]
The stress-strain relationship of the material is simulated by using the yield function f of the equation (1) identified as described above.
[0143]
5 and 6 further show the results of simulating the stress-strain relationship according to the present embodiment for each of the low strain region and the high strain region in solid line graphs. As is apparent from these two graphs, according to the present embodiment, the predicted value by the simulation matches the experimental value with high accuracy not only in the low strain region but also in the high strain region. The Bauschinger effect is accurately simulated even in the high strain region.
[0144]
Further, in the present embodiment, the amount of springback at the time of press forming the material is predicted based on the identified yield function f and the residual stress-strain relationship predicted by the forming simulation. Therefore, according to the present embodiment, since the springback amount is predicted using the yield function f that accurately expresses the Bauschinger effect even in a high strain region, the prediction accuracy is improved.
[0145]
As is clear from the above description, in the present embodiment, step S1 in FIG. 2 constitutes the “experimental value acquisition step”, step S2 constitutes the “coefficient function identification step” and the “relationship identification step”, Steps S3 and S4 cooperate to constitute a “residual stress-strain relationship prediction step”, and steps S5 and S6 together constitute a “spring back amount prediction step”. The equation (1) is a “yield function”, the equation (6) is a “coefficient function” of a coefficient C which is one of the two coefficients of the yield function f, and the equations (7) to (7) 9) is the “coefficient function” of the other coefficient H.
[0146]
In addition, according to the present embodiment, the basic form of the yield function f capable of continuously simulating the stress-strain relationship is predetermined, and a plurality of discrete experimental values are mainly Since it is used to determine various coefficients and constants in the yield function f, the stress-strain relationship can be simulated continuously and with high accuracy from a relatively small number of experimental values.
[0147]
In addition, a program for numerical analysis of the behavior of an elasto-plastic body by a finite element method is generally expressed in increments of tensor quantities. Therefore, in order to incorporate the material hardening model into a finite element method program used in press molding or the like, it is desirable that the material hardening model is also expressed by an increment of the tensor amount. On the other hand, in this embodiment, the composite hardening model, which is the material hardening model, is expressed in increments of the tensor amount, so that the material hardening model can be easily incorporated into the finite element method program. is there.
[0148]
In addition, in this embodiment, in the present embodiment, the “coefficient function identification step” is automatically executed by a computer without substantially requiring human intervention. On the other hand, it can be performed by using a computer as an auxiliary.
[0149]
As mentioned above, although one Embodiment of this invention was described in detail based on drawing, this is an illustration and the aspect described in the above-mentioned section of [the subject which invention is going to solve, a problem-solving means, and the effect of invention] In addition, the present invention can be implemented in various modifications and improvements based on the knowledge of those skilled in the art.
[Brief description of the drawings]
FIG. 1 is a system diagram showing a springback amount prediction apparatus suitable for implementing a springback amount prediction method including a stress-strain relationship simulation method according to an embodiment of the present invention.
FIG. 2 is a process diagram showing the springback amount prediction method.
FIG. 3 is a block diagram for explaining a relationship between three programs executed by the computer 12 of FIG. 1;
FIG. 4 shows the relationship between stress σ and strain ε during tension-compression, and the initial yield stress Y in the composite hardening model used in the embodiment.₀5 is a graph showing the relationship between the back stress α and the isotropic hardening amount R.
FIG. 5 is a graph for comparing the accuracy with which the stress-strain relationship of a material is simulated in a low strain region with reference to experimental values between one prior art proposed by Chaboche-Rousselier and the above embodiment. is there.
FIG. 6 is a graph for comparing the accuracy with which the stress-strain relationship of a material is simulated in a high strain region with reference to experimental values between one prior art proposed by Chaboche-Rousselier and the above embodiment. is there.
FIG. 7 shows the first back stress component α in the embodiment.₁Are analyzed separately for the extremely low strain region (1) and the non-very low strain region (2), and the second back stress component α₂Is a graph showing that the characteristics are analyzed separately in a low strain region (3), a medium strain region (4), and a high strain region (5).
FIG. 8 shows that the proportional coefficient y is equivalent plastic strain ε in the embodiment.^P _eqIt is a graph which shows a mode that it changes according to.
FIG. 9 is a flowchart showing a plastic constitutive equation identification program in FIG. 3;
[Explanation of symbols]
12 Computer
16 External storage device
30 CD-ROM
32 FD

Claims (6)

The stress-strain relationship of a material is approximated by a composite hardening model that combines an isotropic hardening model and a kinematic hardening model. The plastic strain of the material is used as a variable, and the back is defined by the back stress function of the plastic strain. By using the yield function that defines the yield surface of the material by using the stress and the amount of isotropic hardening as variables, the stress-strain relationship of the material can be expressed continuously and the Bauschinger effect of the material can be expressed. A method of simulating
A vector representing a first back stress component expressed by an equation including a coefficient C that is a component that depends on nonlinearity of the kinematic hardening of the back stress and that represents the convergence speed of the non-linear kinematic hardening. , A component that depends on the linearity of kinematic hardening of the back stress and a vector that represents a second back stress component represented by a formula that includes a coefficient H representing the magnitude of the linear kinematic hardening. ,And,
(a) the coefficient C is defined by a coefficient function having a tendency to decrease according to the plastic strain;
(b) The coefficient H is defined by a coefficient function having a tendency to increase according to the plastic strain.
And at least one of
The material is loaded in one direction of the tensile direction and the compression direction and plastically deformed and then unloaded, and further, the material is plastically deformed by loading in the opposite direction, and during that time, the experimental value of the stress-strain relationship of the material Experimental value acquisition step of acquiring a plurality of discrete values;
A coefficient function identifying step of identifying the at least one coefficient function of the coefficient C and the coefficient H based on a plurality of experimental values acquired in the experimental value acquiring step. Relationship simulation method. The stress-strain relationship of a material is approximated by a composite hardening model that combines an isotropic hardening model and a kinematic hardening model. The plastic strain of the material is used as a variable, and the back is defined by the back stress function of the plastic strain. By using the yield function that defines the yield surface of the material by using the stress and the amount of isotropic hardening as variables, the stress-strain relationship of the material can be expressed continuously and the Bauschinger effect of the material can be expressed. A method of simulating
The yield function f is

However,
S _i (i = x, y, z): Deviation stress relating to vertical stress σ _i acting on a plane perpendicular to the coordinate axis corresponding to index i S _ij (i, j = x, y, z): corresponding to index i Deviation stress α _i (i = x, y, z) acting on the shear stress σ _ij acting in the direction of the coordinate axis corresponding to the index j on a plane perpendicular to the coordinate axis to be applied to the plane perpendicular to the coordinate axis corresponding to the index i Back stress α _ij (i, j = x, y, z): Back stress acting in the direction of the coordinate axis corresponding to the index j on a plane perpendicular to the coordinate axis corresponding to the index i Y ₀ : Plastic deformation of the material Initial yield stress, which is the yield stress when starting for the first time, R: expressed by the formula of the amount of isotropic hardening of the material,
A vector {α} representing the back stress α is

However,
{Α ₁ }: A vector representing a first back stress component which is a component depending on nonlinearity of kinematic hardening among back stresses {α ₂ }: a _second component which is a component depending on linearity of kinematic hardening among back stresses A vector {dα ₁ } that is expressed by a vector representing a back stress component and represents an increment dα ₁ of the first back stress component α ₁ is

However,
C: coefficient representing the speed of convergence of nonlinear kinematic hardening 2a / 3: convergence value of non-linear kinematic hardening {dε ^P }: vector representing plastic strain increment dε ^P _i (i = x, y, z) dε ^P _eq : A vector {dα ₂ } expressed by the equivalent plastic strain increment and representing the increment dα ₂ of the second back stress component α ₂ is

However,
H: It is expressed by an expression of a coefficient representing the magnitude of linear kinematic hardening,
(A) The coefficient C is

However,
ε ^P _eq : Equivalent plastic strain C ′, ε ₀ : Expressed by a constant equation, (B) The coefficient H is

However,
ε ₁ : boundary value ε ₂ : boundary value (> ε ₁ )
H ₁ , p, q, H ₂ : including at least one of being expressed by a constant expression, and
The material is loaded in one direction of the tensile direction and the compression direction and plastically deformed and then unloaded, and further, the material is plastically deformed by loading in the opposite direction, and during that time, the experimental value of the stress-strain relationship of the material Experimental value acquisition step of acquiring a plurality of discrete values;
A coefficient function identifying step of identifying the at least one coefficient function of the coefficient C and the coefficient H based on a plurality of experimental values acquired in the experimental value acquiring step. Relationship simulation method. Furthermore, the relationship between the amount of isotropic hardening and the plastic strain is such that the amount of isotropic hardening increases according to the plastic strain regardless of the magnitude of the plastic strain, and the amount of isotropic hardening and the plastic strain. stress according to claim 1 or 2 relationship includes a relationship identification step of identifying, based on the plurality of experimental values obtained - strain relationship simulating method. The increment dR of the isotropic hardening amount R is

However,
b: The stress-strain relationship simulation method according to claim 3 , which is expressed by a constant equation representing the magnitude of isotropic hardening. A method for predicting a springback amount, which is an amount by which a material elastically recovers after the end of plastic working by using the yield function according to any one of claims 1 to 4 ,
A specimen having substantially the same material as the material is loaded in one direction in the tensile direction and the compression direction and plastically deformed, then unloaded, and further loaded in the opposite direction to be plastically deformed. An experimental value acquisition step of acquiring the experimental value of the stress-strain relationship of the test piece as a plurality of discrete values;
The coefficient of the back stress in the yield function is a variable defined by a coefficient function having the plastic strain as a variable, and the coefficient function is identified based on a plurality of acquired experimental values. Process,
A residual stress-strain relationship prediction step for predicting a stress-strain relationship remaining in the material at the end of the plastic working;
A springback amount prediction step for predicting the springback amount based on the predicted residual stress-strain relationship and the yield function including the coefficient defined by the identified coefficient function. To predict the amount of springback to be performed. The recording medium which recorded the program run by a computer in order to implement the springback amount prediction method of Claim 5 , so that computer reading was possible.

Images