JP2004042098A - 成形シミュレーション解析方法 - Google Patents

Abstract

【課題】薄板であるか厚板であるかを問わず、板材のプレス成形による変形をコンピュータによるシミュレーションによって高精度で解析する。
【解決手段】板材のプレス成形による変形をコンピュータによるシミュレーションによって解析するために、板材がそれの面方向に並んだ複数のシェル要素の集合体としてコンピュータ上でモデル化されたシェル要素モデルを用い、各シェル要素の板厚方向応力を考慮した上で各シェル要素の面方向の応力である面内応力を解析する。
【選択図】　図２７

Description

【０００１】
【発明の属する技術分野】
本発明は、板材のプレス成形による変形をコンピュータによるシミュレーションによって解析する技術に関するものである。
【０００２】
【従来の技術】
金型を利用して板材を塑性加工によって成形して製品を製造する分野が存在する。その一例は、自動車を製造する分野である。
【０００３】
板材の成形においては、必要に応じて金型の修正が、製品の成形性、製品の生産性、製品の形状精度、製品の面品質（面精度）等の各要素を改善するために行われる。金型修正に必要な工数および時間の削減が強く要望されており、例えば、成形性および生産性の改善を目的とした金型修正については、近年、工数および時間の削減が着実に達成されつつある。
【０００４】
しかし、形状精度および面品質の改善を目的とした金型修正については、工数および時間の削減がそれほど達成されておらず、特に、形状精度が不良であるという問題が金型修正のための工数の低減を阻害すると同時に生産準備リードタイムの短縮を阻害している。
【０００５】
製品の形状精度は、製品の材料に依存する。例えば、自動車を製造する分野においては、最近、自動車の軽量化および衝突安全性向上のため、製品の材料として高張力鋼板およびアルミニウム合金を使用することが普及しつつある。
【０００６】
一方、製品の形状精度が不良となる原因として、成形後の製品のスプリングバックがある。スプリングバックという現象は、プレス成形等、塑性加工において、材料が負荷された後に除荷されると、材料が弾性回復してその形状が変化する現象である。このスプリングバック量は、製品形状のみならず材料の種類にも依存し、例えば、高張力鋼板およびアルミニウム合金である場合において従来一般に使用された材料である場合より大きいという傾向がある。
【０００７】
形状精度の改善を目的とした金型修正のための工数および時間の削減を図るため、次のような手法が既に提案されている。これは、板材の成形（プレス成形）によって板材が変形する様子を数値シミュレーションによって解析し、その解析結果を金型設計の段階において利用することにより、製品の形状精度不良を事前に予測し、それにより、製品の形状精度が許容範囲内に収まるように金型をコンピュータ上で修正するというものである。
【０００８】
この手法を採用するにしても、実際に製作された金型を使って成形を行って製品の形状精度を測定した後にその形状精度が改善されるようにその金型を機械加工等によって修正する工程が完全に省略されるとは限らない。とはいえ、その手法を採用すれば、金型を仮想的に修正することが可能となり、一方、その修正のための工数および時間は、金型を実際に修正するための工数および時間に比べて容易に削減可能である。したがって、その手法を採用すれば、金型の設計および修正のための工数および時間を全体として削減することが容易となる。
【０００９】
この手法を具体化するため、次のようなスプリングバック量予測方法が既に提案されている。それは、形状精度不良の予測精度を向上させるため、材料モデリングと計算技術との両面から検討を行い、材料の除荷時における材料特性の非線形性をモデル化して数値解析計算としてのＦＥＭ計算に導入したものである。この種のスプリングバック量予測方法は、成形シミュレーション工程とスプリングバック解析工程とを含んでおり、その一従来例が特開２０００−３１２９３３号公報に開示されている。
【００１０】
本発明者らは、この種のスプリングバック量予測方法を具体化して実験を行った結果、薄板鋼板については、金型設計において利用可能な精度でスプリングバック量を予測できることを確認した。
【００１１】
本発明者らは、さらに、同じスプリングバック量予測方法により、厚板についてスプリングバック量を予測する実験を行った。厚板の成形によって製造される製品として、自動車を例にとれば、車体のフレームやサスペンションのアームという機能部品が挙げられる。
【００１２】
【発明が解決しようとする課題】
しかし、本発明者らは、上記の実験により、板材のプレス成形による変形のシミュレーションを、その板材が厚板である場合には薄板である場合ほどには高い精度で行うことが困難であることに気が付いた。
【００１３】
【課題を解決するための手段および発明の効果】
本発明によって下記の各態様が得られる。各態様は、項に区分し、各項に番号を付し、必要に応じて他の項の番号を引用する形式で記載する。これは、本明細書に記載の技術的特徴のいくつかおよびそれらの組合せのいくつかの理解を容易にするためであり、本明細書に記載の技術的特徴やそれらの組合せが以下の態様に限定されると解釈されるべきではない。
【００１４】
（１）　曲げ変形させられるべき板材の両面の一方である曲げ外面をダイスとパッドとのうちの少なくともダイスによって支持しつつ他方の面である曲げ内面にポンチを加圧状態で接触させることによって前記板材をプレス成形する際のその板材の変形をコンピュータによるシミュレーションによって解析する方法であって、
前記板材がそれの面方向に並んだ複数のシェル要素の集合体として前記コンピュータ上でモデル化されたシェル要素モデルを用い、前記各シェル要素の板厚方向応力を考慮した上で各シェル要素の面方向の応力である面内応力を解析する成形シミュレーション工程を含む成形シミュレーション解析方法。
【００１５】
板材のプレス成形による変形をシミュレーションにより予測する技術が既に存在する。
【００１６】
この成形シミュレーションは、板材を、それの面方向と板厚方向との双方に沿って複数のソリッド要素に仮想的に分割してモデル化し、そのソリッド要素モデルのもとにＦＥＭ等の数値解析を行う手法と、板材を、それの面方向に複数のシェル要素に仮想的に分割してモデル化し、そのシェル要素モデルのもとにＦＥＭ等の数値解析を行う手法とに分類される。
【００１７】
一般に、先の手法は、板材変形の予測精度を容易に向上させることが可能である点で後の手法より優れ、一方、後の手法は、数値解析に必要な計算時間を容易に短縮可能な点で先の手法より優れている。したがって、予測精度と計算時間との両方において実用的に優れた手法の開発が要望されている。
【００１８】
後の手法を採用する場合、従来のシェル要素は、ライスナ−ミンドリン（Ｒｅｉｓｓｎｅｒ−Ｍｉｎｄｌｉｎ）の厚板理論に基づく下記の解析条件（ａ）ないし（ｄ）のもとに設計される。
【００１９】
以下、その解析条件を説明するが、その説明の便宜上、板材の中立面の延びる方向すなわち曲げ方向をｘ方向、板厚方向をｚ方向といい、それら両方向に直角な方向をｙ方向（断面方向ともいう）という。それらのうちｘ方向とｙ方向とはいずれも、板材にとっての面内方向であり、残りのｚ方向は、板材の面と交差するという意味において面外方向ということが可能である。
【００２０】
（ａ）シェル要素について平面応力状態を想定し、よって、シェル要素内における板厚方向の垂直応力σｚを無視し、０とみなす。
【００２１】
（ｂ）各シェル要素について平面保持性、すなわち、曲げ変形前にシェル要素軸線に垂直であった断面が曲げ変更後も平面を保持する。したがって、シェル要素の断面回転角θｙ（ｙ軸回りの回転角）は板厚方向の位置すなわちｚ座標値に依存しない。
【００２２】
ここに、断面回転角θｙとは、シェル要素の中立面の法線に対するそのシェル要素の辺の角度を意味する。より正確には、曲げ変形前にシェル要素の長手方向軸線に垂直であった平面断面を変形前断面、曲げ変更後にその変形前断面が移動した断面を変形後断面とした場合に、その変形後断面上の任意の点上の接線と変形前断面との成す角度を意味する。
【００２３】
（ｃ）断面回転角θｙがｚ座標値に依存しないから、シェル要素内における曲げ方向の垂直ひずみεｘは板厚方向の位置に応じて線形に分布し、また、板厚方向のせん断ひずみγｘｚは０でない一定値である。
【００２４】
（ｄ）シェル要素の平面保持性によりそれの中立面は常に板厚中央面と一致する。
【００２５】
この従来のシェル要素を用いた成形シミュレーション解析は、薄板成形のように、板材の曲げの厳しさが比較的弱い場合には、上述の解析条件の採用が妥当であるため、変形を精度よく計算できる。
【００２６】
しかし、本発明者らの研究により、この従来の成形シミュレーションでは厚板の変形を精度よく予測するのは困難であることに気が付いた。厚板の場合には、上述の解析条件を採用することが妥当ではなく、各シェル要素の応力分布（垂直応力σとせん断応力τとによって記述される）および変形状態（垂直ひずみεとせん断ひずみγとによって記述される）をより高精度に再現する新たな解析モデルおよび解析技術の確立が必要である。
【００２７】
そのためには、板厚方向応力σｚを考慮する場合と無視する場合とで曲げ方向応力σｘ（断面方向応力σｙも同様）の計算値が互いに一致しないという事実に着目し、板厚方向応力σｚを考慮して曲げ方向応力σｘ（断面方向応力σｙも同様）を計算することが重要である。
【００２８】
この知見に基づき、本項に係る方法によれば、板材がそれの面方向に並んだ複数のシェル要素の集合体としてコンピュータ上でモデル化されたシェル要素モデルを用い、各シェル要素の板厚方向応力を考慮した上で各シェル要素の面方向の応力である面内応力が解析され、その面内応力の解析結果を用いて板材の変形が解析される。
【００２９】
このように、この方法によれば、板材の変形がシェル要素モデルを用いて解析されるから、ソリッド要素モデルを用いる場合に比較し、板材の変形をシミュレートするための計算時間を短縮することが容易となる。
【００３０】
しかも、この方法によれば、前述の従来のシェル要素モデルを用いて板材の変形を解析する場合とは異なり、その解析が板厚方向応力を考慮して行われるため、その解析の精度を、ソリッド要素モデルを用いる場合に近づけることが容易となる。
【００３１】
このように、この方法によれば、計算時間の短縮と解析精度の向上とを両立させることが容易となる。
【００３２】
この方法の用途の一例は、成形シミュレーションの解析結果に基づいて板材のスプリングバックを解析することであるが、この方法は他の用途に適用することが可能である。
【００３３】
本項および下記の各項に係る方法は、板材が厚板である場合に有効であるのはもちろんであるが、厚板の成形シミュレーションのみならず薄板の成形シミュレーションにも適用することが可能である。
【００３４】
本項および下記の各項において「応力」という用語は、垂直応力とせん断応力との少なくとも一方を包含し、「面内応力」という用語は、曲げ方向応力と断面方向応力との少なくとも一方を包含し、「ひずみ」という用語は、垂直ひずみとせん断ひずみとの少なくとも一方を包含する。
【００３５】
本項における「成形シミュレーション工程」は、板厚方向応力を考慮して面内応力を解析した後、面内応力と板厚方向応力とのうち、面内応力は用いるが板厚方向応力は用いないで板材の変形を解析する態様で実施したり、それら面内応力と板厚方向応力との双方を用いて板材の変形を解析する態様で実施することが可能である。前者の態様は、後者の態様より単純なアルゴリズムで板材の変形を解析することが容易である。
【００３６】
（２）　前記成形シミュレーション工程が、
前記各シェル要素につき、前記板厚方向応力が０である平面応力状態を想定して前記面内応力の暫定値を解析する暫定値解析工程と、
前記各シェル要素につき、前記解析された暫定値に基づき、かつ、前記板厚方向における力のつり合いを考慮することにより、前記面内応力の最終値を解析する最終値解析工程と
を含む（１）項に記載の成形シミュレーション解析方法。
【００３７】
シェル要素モデルを利用する場合には、板厚方向応力σｚを直接に計算することは困難であるが、例えば部分円筒を成す微小要素を想定するなどして、曲げ変形させられるべき板材における板厚方向（半径方向）における力のつり合いを式で表現すれば、板厚方向応力σｚと曲げ方向応力σｘ（断面方向応力σｙも同様）との関係が誘導できる。したがって、その関係を利用すれば、シェル要素モデルを用いて計算された曲げ方向応力σｘ（断面方向応力σｙも同様）から板厚方向応力σｚを計算することが可能となる。
【００３８】
このような知見に基づき、本項に係る方法によれば、各シェル要素につき、板厚方向応力が０である平面応力状態を想定して面内応力の暫定値が解析される。さらに、各シェル要素につき、その解析された暫定値に基づき、かつ、板厚方向における力のつり合いを考慮することにより、面内応力の最終値が解析される。
【００３９】
したがって、この方法によれば、板厚方向応力を直接に計算することは困難であるシェル要素モデルを用いるにもかからわず、そのシェル要素モデルを用いて計算された面内応力を有効に用いることにより、板厚方向応力を考慮した面内応力を計算することが可能となる。
【００４０】
（３）　前記最終値解析工程が、前記面内応力の、前記各シェル要素の断面上における総和である断面力が、その断面に隣接したシェル要素についての断面力と実質的に一致するように前記最終値を決定する最終値決定工程を含む（２）項に記載の成形シミュレーション解析方法。
【００４１】
前記（２）項に係る方法は、各シェル要素につき、板厚方向における力のつり合いを考慮することにより、面内応力の最終値を、それの暫定値とは異なる値として決定する態様で実施することが可能である。
【００４２】
各シェル要素につき、面内応力の暫定値、すなわち、平面応力状態において板厚方向応力を無視して計算された面内応力とは異なる値を結果的に有するように面内応力の最終値を決定することは、各シェル要素の断面上における面内応力の総和である断面力が、その断面に隣接したシェル要素についての断面力と異なってしまう可能性がある。このような事態は、互いに隣接したシェル要素により共有される節点における面内方向のつり合いに反する。
【００４３】
このような知見に基づき、本項に係る方法によれば、面内応力の、各シェル要素の断面上における総和である断面力が、その断面に隣接したシェル要素についての断面力と実質的に一致するように、面内応力の最終値が決定される。
【００４４】
したがって、この方法によれば、平面応力状態で計算された面内応力が板厚方向応力を考慮して変更されるにもかかわらず、互いに隣接したシェル要素間において断面力の連続性を保持することが可能となる。
【００４５】
（４）　前記最終値決定工程が、前記各シェル要素の断面上において、前記板厚方向応力を考慮した面内応力の総和である断面力が前記平面応力状態における面内応力の総和である断面力に実質的に維持されるように前記各シェル要素の中立面の位置を決定し、その決定された中立面位置に基づき、前記最終値を決定する工程を含む（３）項に記載の成形シミュレーション解析方法。
【００４６】
この方法によれば、各シェル要素の断面に作用する面内応力の総和である断面力が実質的に維持されるように中立面の位置が決定され、その結果、互いに隣接したシェル要素間において断面力の連続性を保持することが可能となる。
【００４７】
（５）　前記最終値解析工程が、前記解析された暫定値と、前記プレス成形時に前記板材が前記曲げ内面において前記ポンチから受ける接触圧とに基づき、前記最終値を解析する工程を含む（２）ないし（４）項のいずれかに記載の成形シミュレーション解析方法。
【００４８】
本発明者らは、板材のプレス成形による変形を精度よくシミュレートするために、成形時にポンチが板材の曲げ内面に接触する接触力を考慮して面内応力を解析することが重要であることに気が付いた。
【００４９】
この知見に基づき、本項に係る方法によれば、面内応力の暫定値と、プレス成形時に板材が曲げ内面においてポンチから受ける接触圧とに基づき、面内応力の最終値が解析される。
【００５０】
（６）　前記成形シミュレーション工程が、
前記各シェル要素につき、前記板厚方向応力が０である平面応力状態を想定して前記面内応力の暫定値を解析する暫定値解析工程と、
前記各シェル要素につき、前記解析された暫定値に基づき、前記板厚方向応力の前記板厚方向における分布を近似的に定義する応力近似関数を特定し、その特定された応力近似関数を用いて前記面内応力の最終値を解析する最終値解析工程と
を含む（１）項に記載の成形シミュレーション解析方法。
【００５１】
この方法は、前記（２）項に係る方法に対し、それの一実施態様を構成する関係にあると考えることが可能である。具体的には、本項に係る方法が実施されることは、板厚方向応力の板厚方向における分布を近似的に定義する応力近似関数を特定し、その特定された応力近似関数を用いて面内応力の最終値を解析する態様で、前記（２）項に係る方法が実施されることに等しい。
【００５２】
（７）　前記最終値解析工程が、
前記解析された暫定値に基づき、前記応力近似関数を特定する応力近似関数特定工程と、
その特定された応力近似関数を用い、前記板厚方向応力と共に成立する面内応力を仮面内応力として決定する仮面内応力決定工程と、
その決定された仮面内応力に基づき、前記面内応力の、前記各シェル要素の断面上における総和である断面力が、その断面に隣接したシェル要素についての断面力と実質的に一致するように前記最終値を決定する最終値決定工程と
を含む（６）項に記載の成形シミュレーション解析方法。
【００５３】
この方法によれば、応力近似関数を用いることにより、板厚方向応力を考慮しない面内応力が、板厚方向応力を考慮した面内応力すなわち仮面内応力に変更される。さらに、その仮面内応力に基づき、シェル要素間における力のつり合いを考慮することにより、面内応力の最終値が決定される。
【００５４】
したがって、この方法によれば、板厚方向応力を考慮しない面内応力を暫定値とし、板厚方向応力と、シェル要素間における力のつり合いとの双方を考慮することにより、面内応力の最終値が決定される。
【００５５】
（８）　前記補正工程が、前記各シェル要素の断面上において、前記板厚方向応力を考慮した面内応力の総和である断面力が前記平面応力状態における面内応力の総和である断面力に実質的に維持されるように前記各シェル要素の中立面の位置を決定し、その決定された中立面位置に基づき、前記最終値を決定する工程を含む（７）項に記載の成形シミュレーション解析方法。
【００５６】
この方法は、前記（４）項に係る方法が前記（３）項に係る方法に対して有する関係と同じ関係を前記（７）項に係る方法との間において有する。
【００５７】
（９）　前記最終値解析工程が、前記解析された暫定値と、前記プレス成形時に前記板材が前記曲げ内面において前記ポンチから受ける接触圧とに基づき、前記応力近似関数を特定する工程を含む（２）ないし（８）項のいずれかに記載の成形シミュレーション解析方法。
【００５８】
本発明者らは、板厚方向応力の板厚方向における分布を近似的に定義する応力近似関数は、面内応力のみを考慮して特定することは可能であるが、面内応力のみならず、成形時にポンチが板材の曲げ内面に接触する接触力をも考慮して特定すれば、板材のプレス成形による変形の解析精度を容易に向上できることに気が付いた。
【００５９】
このような知見に基づき、本項に係る方法によれば、面内応力の暫定値と、プレス成形時に板材が曲げ内面においてポンチから受ける接触圧とに基づき、板厚方向応力の応力近似関数が特定される。
【００６０】
（１０）　前記応力近似関数が、前記各シェル要素につき、前記板厚方向における力のつり合いを前提として前記分布を定式化する（６）ないし（９）項のいずれかに記載の成形シミュレーション解析方法。
【００６１】
（１１）　前記各シェル要素が複数の節点により定義され、かつ、それら複数の節点が、対応するシェル要素を幾何学的に特定するのに不可欠ではない節点を含まない（１）ないし（１０）項のいずれかに記載の成形シミュレーション解析方法。
【００６２】
この方法によれば、各シェル要素を定義する複数の節点が、同じシェル要素を幾何学的に特定するのに必要最小の節点により構成されることとなる。したがって、この方法によれば、板材全体に設定しなければならない節点の数を容易に削減可能となり、ひいては節点数の削減による計算時間の短縮が容易となる。
【００６３】
（１２）　前記成形シミュレーション工程が、前記プレス成形による曲げの厳しさに関して前記板材と等価な等価板材がそれの面方向と板厚方向との双方に並んだ複数のソリッド要素の集合体として前記コンピュータ上でモデル化されたソリッド要素モデル（またはソリッド要素モデルを用いた計算結果と、モデル実験による測定結果）を用い、前記曲げ変形前に前記等価板材の長手方向軸線に垂直である平面断面上の各点が前記曲げ変形によって回転する角度である断面回転角を前記等価板材の板厚方向位置に依存するように解析し、その断面回転角の解析結果を用い、前記等価板材のひずみを解析するひずみ解析工程を含む（１）ないし（１１）項のいずれかに記載の成形シミュレーション解析方法。
【００６４】
本発明者らは、種々の研究の結果、厚板成形においては、薄板成形とは異なり、各シェル要素について平面保持性が成立しないから、断面回転角θｙをｘ座標値のみならずｚ座標値にも依存させることが重要であることに気が付いた。
【００６５】
本発明者らは、さらに、断面回転角θｙをｚ座標値に依存させることにより、ひずみの板厚方向における非線形性分布を考慮することが重要であることにも気が付いた。
【００６６】
本発明者らは、さらにまた、プレス成形による曲げの厳しさ（例えば、板材の曲げ半径Ｒを板厚ｔで割り算した値である曲げ指標Ｒ／ｔによる表現することができる）が同じである限り、曲げ半径Ｒまたは板厚ｔが異なっても、断面回転角およびひずみの板厚方向における分布が共通するという事実にも気が付いた。
【００６７】
したがって、複数の板材について成形シミュレーションを行うことが必要である場合に、上述のソリッド要素モデルを用いた断面回転角の解析（またはモデル実験）を各板材ごとに行うことは不可欠ではない。
【００６８】
それらの知見に基づき、本項に係る方法によれば、プレス成形による曲げの厳しさに関して、成形シミュレーションの実行対象である板材と等価な等価板材がそれの面方向と板厚方向との双方に並んだ複数のソリッド要素の集合体としてコンピュータ上でモデル化されたソリッド要素モデル（またはソリッド要素モデルを用いた計算結果と、モデル実験による測定結果）を用い、曲げ変形前に等価板材の長手方向軸線に垂直である平面断面上の各点が曲げ変形によって回転する角度である断面回転角が等価板材の板厚方向位置に依存するように解析される。さらに、その断面回転角の解析結果を用いて等価板材のひずみが解析される。
【００６９】
本項において「等価板材」は、面内応力を計算すべき板材そのものを排除しない用語として使用される。また、本項において「ひずみ」は、曲げ方向ひずみと板厚方向ひずみとの少なくとも一方を含み、また、各方向のひずみは、垂直ひずみとせん断ひずみとの少なくとも一方を含む。
【００７０】
（１３）　前記曲げの厳しさが、前記曲げ変形後の前記板材の曲げ半径が小さいほど強くなり、かつ、前記板材の板厚が厚いほど強くなるように定義される（１２）項に記載の成形シミュレーション解析方法。
【００７１】
（１４）　前記曲げの厳しさが、前記曲げ変形後の前記板材の曲げ半径をその板材の板厚で割り算した曲げ指標により表現される（１３）項に記載の成形シミュレーション解析方法。
【００７２】
（１５）　当該方法が、前記曲げ指標が３以下である板材について実施される（１４）項に記載の成形シミュレーション解析方法。
【００７３】
（１６）　前記ひずみ解析工程が、前記等価板材につき、前記断面回転角の解析結果に基づき、前記ひずみの前記板厚方向における分布を近似的に定義するひずみ近似関数を特定するひずみ近似関数特定工程を含み、
前記成形シミュレーション工程が、前記各シェル要素につき、前記シェル要素モデルを前記特定されたひずみ近似関数と共に用い、かつ、前記板厚方向応力を考慮した上で前記面内応力を解析する面内応力解析工程を含む（１２）項に記載の成形シミュレーション解析方法。
【００７４】
この方法によれば、前記ソリッド要素モデル（またはソリッド要素モデルを用いた計算結果と、モデル実験による測定結果）を用いて曲げの厳しさに関連付けて計算された断面回転角に基づいてひずみ近似関数が特定され、その特定されたひずみ近似関数と共に前記シェル要素モデルを用いることにより、面内応力が解析される。
【００７５】
したがって、この方法によれば、面内応力を解析するためにシェル要素モデルを用いるにもかかわらず、ひずみの板厚方向における非線形分布を考慮して面内応力を解析することが容易となる。
【００７６】
（１７）　さらに、前記板材をプレス成形によって塑性加工した後にその板材を除荷することに伴ってその板材が弾性回復する現象であるスプリングバックをコンピュータによって解析するスプリングバック解析工程を含み、
前記成形シミュレーション工程が、前記除荷時に前記板材に残留する残留応力と、前記除荷時に前記板材に残留する残留ひずみとを解析するものであり、
前記スプリングバック解析工程が、それら残留応力および残留ひずみの各解析値を用い、前記スプリングバックを解析するものである（１）ないし（１６）項のいずれかに記載の成形シミュレーション解析方法。
【００７７】
この方法によれば、板厚方向応力を考慮して高精度で解析された面内応力を利用することにより、板材のスプリングバックを精度よく解析することが可能となる。
【００７８】
（１８）　さらに、
前記板材と材質が実質的に同じ試験片を引張方向と圧縮方向との一方向に負荷して塑性変形させた後に除荷し、さらに、逆方向に負荷して塑性変形させるとともに、その間にその試験片の応力−ひずみ関係について実験値を取得することを、試験片の除荷開始時における塑性ひずみである予ひずみを変化させるごとに繰り返す実験値取得工程と、
その取得された実験値により表される前記応力−ひずみ関係を線形の材料硬化モデルで近似することにより、前記スプリングバックを解析するのに必要な材料特性値を前記予ひずみに関連付けて取得する材料特性値取得工程と
を含み、かつ、
前記スプリングバック解析工程が、前記残留ひずみの解析値を前記予ひずみの今回値として見た場合に、その予ひずみの今回値に関連付けて取得された材料特性値を用いるとともに、前記残留ひずみおよび残留応力の各解析値を用い、前記線形の材料硬化モデルに従って前記スプリングバックを解析するものである（１７）項に記載の成形シミュレーション解析方法。
【００７９】
（１９）　（１）ないし（１８）項のいずれかに記載の成形シミュレーション解析方法を実施するためにコンピュータにより実行されるプログラム。
【００８０】
このプログラムがコンピュータにより実行されれば、前記（１）ないし（１８）項のいずれかに係る方法と同じ作用効果が実現され得る。
【００８１】
本項における「プログラム」は、それの機能を果たすためにコンピュータにより実行される指令の組合せのみならず、各指令により処理されるファイルやデータをも含むように解釈することが可能である。
【００８２】
また、本項における「プログラム」は、各項に係る機能を果たすためにあえて作成された専用のプログラムとして、または専用のサブプログラム（モジュール等）の組合せとして作成することは不可欠ではなく、専用のサブプログラムと汎用のサブプログラムとの組合せとして作成することが可能であり、場合によっては、専用のサブプログラムを伴わない複数の汎用のサブプログラムの組合せとして作成する場合もあり得る。
【００８３】
（２０）　（１９）項に記載のプログラムをコンピュータ読み取り可能に記録した記録媒体。
【００８４】
この記録媒体に記録されているプログラムがコンピュータにより実行されれば、前記（１）ないし（１８）項のいずれかに係る方法と同じ作用効果が実現され得る。
【００８５】
本項における「記録媒体」は種々の形式を採用可能であり、例えば、フレキシブル・ディスク等の磁気記録媒体、ＣＤ、ＣＤ−ＲＯＭ等の光記録媒体、ＭＯ等の光磁気記録媒体、ＲＯＭ等のアンリムーバブル・ストレージ等の少なくとも１つを採用可能である。
【００８６】
【発明の実施の形態】
以下、本発明のさらに具体的な実施の形態の一つを図面に基づいて詳細に説明する。
【００８７】
本発明の一実施形態はスプリングバック解析方法であり、図１には、そのスプリングバック解析方法を実施するのに好適なスプリングバック解析装置（以下、単に「解析装置」という）が示されている。
【００８８】
この解析装置は、入力装置１０とコンピュータ１２と出力装置１４と外部記憶装置１６とを備えている。入力装置１０はマウス，キーボード等を含むように構成される。コンピュータ１２はＣＰＵ等、プロセッサ２０と、ＲＯＭ，ＲＡＭ，ハードディスク等、メモリ２２と、それらプロセッサ２０とメモリ２２とを互いに接続するバス２４とを含むように構成される。出力装置１４はディスプレイ，プリンタ，プロッタ等として構成される。
【００８９】
外部記憶装置１６は、ＣＤ−ＲＯＭ３０，書き込み可能なＭＯ（光磁気ディスク）３２等、記録媒体が装填可能となっていて、装填状態においては、記録媒体に対するデータの読み取りおよび書き込みが必要に応じて行われる。
【００９０】
本実施形態においては、スプリングバックの解析に必要なプログラムがＣＤ−ＲＯＭ３０に記憶され、また、その解析に必要なデータがＭＯ３２に記憶されている。スプリングバックの解析に必要なプログラムには、材料特性値取得プログラムと、成形シミュレーション・プログラムと、材料特性値決定プログラムと、バックストレス計算プログラムと、スプリングバック解析プログラムとがある。
【００９１】
スプリングバックの解析時には、それらＣＤ−ＲＯＭ３０およびＭＯ３２から必要なプログラムおよびデータが読み出されてコンピュータ１２のＲＡＭまたはハードディスクに転送され、その後、プロセッサ２０によりそのプログラムが実行される。
【００９２】
図２には、本実施形態であるスプリングバック解析方法が工程図で表されている。
【００９３】
概略的に説明すれば、スプリングバック解析方法は、板材としての厚板鋼板をプレス成形した後にその板材に生じるスプリングバックを解析するものである。そのプレス成形は、塑性ひずみが０．１５，０．２０というように、比較的高い塑性ひずみが生じる部分が存在するように行われる。
【００９４】
なお付言すれば、以下、「応力」という用語は、特に断りがない限り、垂直応力とせん断応力との双方を意味し、同様に、「ひずみ」という用語も、垂直ひずみとせん断ひずみとの双方を意味する。
【００９５】
図３には、プレス成形機１００が例示されている。このプレス成形機１００は、ワークＷとしての厚板をしわ押さえを利用せずにＵ字状に曲げ加工することが可能なものである。このプレス成形機１００においては、上ラム１０２と下ラム１０４とがすきまを隔てて互いに対向させられ、それら上ラム１０２と下ラム１０４との間において、ポンチ１１０とダイス１１２とパッド１１４とがワークＷを挟むように配置される。
【００９６】
ポンチ１１０は、ダイス１１２に接近する向きの力を上ラム１０２から受けて、予め設定された速度でワークＷに押し付けられる。ダイス１１２はステージ１２０に固定されている。パッド１１４は、油圧シリンダ１２２を介して下ラム１０４に支持されていて、ポンチ１１０に接近する向きの力を油圧シリンダ１２２から受けるようにされている。ダイス１１２はダイセット１２４によってステージ１２０に取り付けられる。
【００９７】
本実施形態のスプリングバック解析方法においては、プレス成形されるべき材料と材質が実質的に同じ試験片を用いることにより、その材料の応力−ひずみ関係の実験値が複数の離散値として取得され、その取得された複数の実験値に基づき、その材料のスプリングバックを解析するのに必要な材料特性値が取得される。この材料特性値は、材料のバウシンガ効果が精度よく表現されるようにその材料の応力−ひずみ関係を定義するのに必要な情報であるということもできる。
【００９８】
また、このスプリングバック解析方法においては、プレス成形直後に材料に残存する応力−ひずみ関係を解析するために、有限要素法を用いることにより、材料のプレス成形による変形をシミュレートする成形シミュレーションが行われる。
【００９９】
このスプリングバック解析方法においては、その後、その有限要素法を用いるとともに、
（１）材料特性値取得プログラムの実行により取得された材料特性値と、
（２）成形シミュレーション・プログラムの実行により計算された残留応力および残留ひずみ、板材の各要素ごとの中立面位置、ならびに板材の各要素ごとの板厚と、
（３）バックストレス計算プログラムの実行により計算されたバックストレスとを用いることにより、材料のスプリングバックが解析される。
【０１００】
さらに具体的に説明すれば、本実施形態においては、図２に示すように、まず、ステップＳ１において、作業者により、プレス成形されるべき材料と材質が実質的に同じ試験片が準備され、その後、図示しない試験機を用いることにより、その試験片に対して一軸応力状態において引張−圧縮試験が行われる。
【０１０１】
この試験は具体的には、試験片を引張方向に負荷して塑性変形させた後に除荷し、さらに、逆方向すなわち圧縮方向に負荷して塑性変形させるというものであり、その間にその試験片の応力−ひずみ関係の実験値が複数の離散値として取得される。
【０１０２】
それら複数の実験値は、試験片の除荷開始時における塑性ひずみである予ひずみを変化させるごとに取得される。図４には、除荷および圧縮時における応力−ひずみ関係を表すグラフが５つの予ひずみについてそれぞれ示されている。
【０１０３】
次に、ステップＳ２において、コンピュータ１２により、前記材料特性値取得プログラムが実行され、それにより、上記取得された複数の実験値に基づいて材料特性値が各予ひずみごとに決定される。材料特性値取得プログラムの内容については、後に詳述する。
【０１０４】
その後、ステップＳ３において、作業者により、前記成形シミュレーションを実行するのに必要なデータ、すなわち、数値解析用データが作成される。数値解析用データは、プレス成形すべき板材を仮想的に複数のシェル要素に分割するためのデータと、プレス成形時に板材に付与される荷重（板材にポンチおよびパッドから負荷される荷重を含む）に関する条件と、複数のシェル要素に関する境界条件とを含んでいる。
【０１０５】
本実施形態においては、各シェル要素を定義する複数の節点が、各シェル要素を幾何学的に特定するのに不可欠ではない節点は含まないようにされている。このようなシェル要素の一例は、四辺形を成し、それの４つの頂点を４つの節点とする四辺形要素である。各シェル要素を上記のように定義することにより、板材に設定される節点の数の増加が回避され、ひいては、節点数の増加に伴う計算時間の増加も回避される。
【０１０６】
続いて、ステップＳ４において、コンピュータ１２により、上記成形シミュレーション・プログラムが実行され、それにより、プレス成形直後に板材に残存する応力−ひずみ関係が解析される。さらに、プレス成形直後における材料の形状も解析される。
【０１０７】
この成形シミュレーション・プログラムは、概略的に説明すれば、材料の応力−ひずみ関係を近似する材料硬化モデルとして等方硬化モデルを使用する形式である。この成形シミュレーション・プログラムは、板材の各シェル要素の番号ｊに関連付けて、プレス成形直後に各シェル要素に存在する塑性ひずみである初期ひずみを残留ひずみとしてメモリ２２に格納する。この成形シミュレーション・プログラムの詳細については、後述する。
【０１０８】
この成形シミュレーション・プログラムは、ＬＳＴＣ社により製造され、「ＬＳ−ＤＹＮＡ３Ｄ」という名称で日本総合研究所により販売されたプログラム、およびＥＳＩ社により製造され、「ＰＡＭ−ＳＴＡＭＰ」という名称で販売されたプログラムを用いるとともに、それらに新たなステップが追加されて構成されている。
【０１０９】
その後、ステップＳ５において、コンピュータ１２により、前記材料特性値決定プログラムが実行され、それにより、上記ステップＳ２において取得された複数の材料特性値を用いることにより、板材の各シェル要素の成形に伴う塑性ひずみに対応する材料特性値が決定される。材料特性値決定プログラムの内容については、後に詳述する。
【０１１０】
続いて、ステップＳ６において、コンピュータ１２により、前記バックストレス計算プログラムが実行され、それにより、板材の各シェル要素ごとにバックストレスαが計算される。バックストレス計算プログラムの内容については、後に詳述する。
【０１１１】
その後、ステップＳ７において、コンピュータ１２により、前記スプリングバック解析プログラムが実行され、それにより、スプリングバックが解析される。このスプリングバック解析プログラムは、前述の説明から明らかなように、
（１）前記材料特性値決定プログラムの実行により板材の各シェル要素ごとに決定された材料特性値と、
（２）前記成形シミュレーション・プログラムの実行により各シェル要素ごとに計算された残留応力および残留ひずみ、中立面位置ｒ_０ならびに板厚ｔ_０と、
（３）前記バックストレス計算プログラムの実行により各シェル要素ごとに計算されたバックストレスと
を用いることにより、板材のスプリングバックを解析するためのプログラムである。
【０１１２】
このスプリングバック解析プログラムの一例は、日本総合研究所により製造されて「ＪＯＨ−ＮＩＫＥ３Ｄ」という名称で販売されたものである。
【０１１３】
以上で、スプリングバック解析方法の一回の実行が終了する。
【０１１４】
前記スプリングバック解析プログラムは、スプリングバックを解析すべき材料の応力−ひずみ関係を種々の方式で定義可能に設計されている。その方式の一つに、応力−ひずみ関係を、平行四辺形を用いた線形の材料硬化モデルで近似する方式がある。図５には、試験片を引張方向に負荷（正負荷）して塑性変形させた後に除荷し、引き続いて逆向きすなわち圧縮方向に負荷（逆負荷）して再度塑性変形させた場合に、塑性ひずみε^Ｐと応力σとが変化する様子が平行四辺形で示されている。
【０１１５】
ここで、図５における各種記号を説明すれば、「Ｙ_０」は、材料が初めての塑性変形を開始するときの降伏応力である初期降伏応力を示している。「Ｙ_１」は、材料が再降伏するときの応力である再降伏応力を示している。「Ｈ」は、塑性硬化係数を示している。塑性硬化係数Ｈは、正負荷時および逆負荷時の塑性域における応力−ひずみ関係を表す直線グラフ（後述の第３直線部５４）の勾配を角度φで表した場合にｔａｎφで表される。「β」は、再降伏応力Ｙ_１が初期降伏応力Ｙ_０に対して低下する程度を規定する材料軟化係数を示している。「ε^Ｐ _０　」は、除荷開始時における塑性ひずみε^Ｐである初期ひずみを示している。
【０１１６】
スプリングバック解析プログラムは、応力−ひずみ関係を近似する材料硬化モデルとして、等方硬化モデルと移動硬化モデルと複合硬化モデルとの中から適当なものを選択できる。
【０１１７】
この選択は、上記材料軟化係数βの値を決定することによって行われる。具体的には、図５に示すように、β＝０のときには、移動硬化モデルが選択され、このとき、再降伏応力Ｙ_１は、
Ｙ_１＝−Ｙ_０＋Ｈ・ε^Ｐ _０
なる式で定義される。また、β＝１のときには、等方硬化モデルが選択され、このとき、再降伏応力Ｙ_１　は、
Ｙ_１＝−Ｙ_０−Ｈ・ε^Ｐ _０
なる式で定義される。また、０＜β＜１のときには、複合硬化モデルが選択され、このとき、再降伏応力Ｙ_１　は、
Ｙ_１＝−Ｙ_０＋Ｈ・（１−２β）ε^Ｐ _０
なる式で定義される。
【０１１８】
図６のグラフは、プレス成形されるべき材料に関するものであるが、試験片に関しても同じグラフが描かれることになる。そして、試験片の応力−ひずみ関係を表すグラフを平行四辺形で近似する場合には、試験片の、正負荷時の弾性域における応力−ひずみ関係を表す第１直線部５０と、除荷時または逆負荷時の弾性域における応力−ひずみ関係を表す第２直線部５２とが互いに平行となるとともに、試験片の、正負荷時の塑性域における応力−ひずみ関係を表す第３直線部５４と、逆負荷時の塑性域における応力−ひずみ関係を表す第４直線部５６とが互いに平行となるように応力−ひずみ関係が近似される。
【０１１９】
一方、応力−ひずみ関係を表すグラフのうちバウシンガ効果が現れる領域は、除荷開始点Ｐ_Ｓ　から再降伏点Ｐ_Ｙ　までの弾性域と、その再降伏点Ｐ_Ｙ　から逆負荷終了点までの塑性域とである。
【０１２０】
したがって、本実施形態においては、それら２つの領域における応力−ひずみ関係を近似的に表す折れ線、すなわち、第２直線部５２と第４直線部５６とが互いに連結されたものの幾何学的特徴が取得されるとともに、その取得された幾何学的特徴に基づき、直接的には試験片、間接的には材料につき、初期降伏応力Ｙ_０と、再降伏応力Ｙ_１と、塑性硬化係数Ｈと、材料軟化係数βとが取得される。
【０１２１】
具体的には、再降伏応力Ｙ_１は、上記折れ線の折れ点における応力σとして取得される。塑性硬化係数Ｈは、第３直線部５４が第４直線部５６と平行であることを利用することにより、第４直線部５６の勾配を角度φで取得してそれのｔａｎ　φを算出することによって取得される。材料軟化係数βは、初期降伏応力Ｙ_０と、再降伏応力Ｙ_１と、塑性硬化係数Ｈと、予ひずみε^Ｐ _ｇとを、前述の式に類似した式、すなわち、
Ｙ_１＝−Ｙ_０＋Ｈ・（１−２β）ε^Ｐ _ｇ
なる式に代入してβについて解くことにより取得される。なお、初期降伏応力Ｙ_０は、正負荷時における実験値に基づいて取得される。
【０１２２】
このスプリングバック解析プログラムは、それら初期降伏応力Ｙ_０，再降伏応力Ｙ_１，塑性硬化係数Ｈおよび材料軟化係数βの他に、材料の弾性係数Ｅ（ヤング率Ｅとせん断弾性係数Ｇとを選択的に表す）を利用することにより、スプリングバックを解析するように設計されている。すなわち、このスプリングバック解析プログラムにおいては、それら初期降伏応力Ｙ_０，再降伏応力Ｙ_１，塑性硬化係数Ｈ，材料軟化係数βおよび弾性係数Ｅがそれぞれ「材料特性値」を構成しているのである。
【０１２３】
ところで、図５のグラフの横軸に取られているのは、弾性ひずみと塑性ひずみとの和である合計ひずみεではなく、塑性ひずみε^Ｐである。したがって、同図のグラフは、例えば、正負荷時の弾性域において、応力σを表す軸に平行に延びている。これに対して、図６のグラフの横軸に取られているのは、合計ひずみεである。したがって、同図のグラフは、例えば、正負荷時の弾性域において、応力σを表す軸に対して傾斜している。よって、それら２つのグラフにおいて、第１ないし第４直線部５０，５２，５４，５６の対応関係はそれら２つの図に示すものとなる。
【０１２４】
したがって、図６における第１直線部５０の勾配を角度γで表した場合にｔａｎγとして求めた値は、材料のヤング率Ｅを表すことになる。一方、第１直線部５０と第２直線部５２とが互いに平行であると仮定されている。したがって、本実施形態においては、第２直線部５２の勾配からヤング率Ｅが取得される。
【０１２５】
その取得されたヤング率Ｅから、予め定められた規則に従い、同じ材料についてのせん断弾性係数Ｇが誘導される。これにより、同じ材料につき、弾性係数Ｅと総称されるヤング率Ｅとせん断弾性係数Ｇとが取得されることになる。
【０１２６】
応力−ひずみ関係は、材料の初期ひずみε^Ｐ _０（試験片の予ひずみε^Ｐ _ｇに相当する）によって変化する。応力−ひずみ関係を近似的に表す平行四辺形の幾何学的特徴が初期ひずみε^Ｐ _０によって変化するのであり、よって、上記材料特性値も初期ひずみε^Ｐ _０によって変化する。したがって、本実施形態においては、前述のように、試験片の予ひずみε^Ｐ _ｇを変化させるごとに引張−圧縮試験が繰り返される。
【０１２７】
図７には、本実施形態で使用される線形の材料硬化モデルが実線グラフで示され、等方硬化モデルが二点鎖線グラフで示され、実際の応力−ひずみ関係が破線グラフで示されている。その実際の応力−ひずみ関係は、正負荷時の弾性域および塑性域においては等方硬化モデルで表される関係（二点鎖線）と一致し、逆負荷時の弾性域および初期の塑性域においては本実施形態の材料硬化モデルで表される関係（実線）と一致する。本実施形態で使用される材料硬化モデルは、線形であるが、予ひずみε^Ｐ _ｇに応じて形状が変化するように定義されている。なお、本実施形態で使用される線形の移動硬化モデルを表す平行四辺形は、同図においては、それの２つの鋭角を含む一対の部分が省略されている。
【０１２８】
本実施形態においては、材料硬化モデルとして、バウシンガ効果を表現し得る移動硬化モデルまたは複合硬化モデルを使用することができるとともに、材料の各シェル要素について予想される初期ひずみε^Ｐ _０に応じて形状が異なる材料硬化モデルが使用される。ただし、各シェル要素について予想される初期ひずみε^Ｐ _０が異なっても、材料硬化モデルが平行四辺形で表現されることは維持される。
【０１２９】
したがって、図７に示すように、本実施形態で使用される線形の材料硬化モデルは、実際の応力−ひずみ関係を表すグラフと全体において一致するわけではないが、バウシンガ効果が現れる領域、すなわち、除荷後に弾性変形する領域と塑性変形する領域とについては、十分に高い精度で一致する。
【０１３０】
このような理由から、本実施形態においては、材料の応力−ひずみ関係が、平行四辺形を用いた線形の材料硬化モデルで近似させられているのであり、よって、本実施形態によれば、低ひずみ領域のみならず高ひずみ領域においてもバウシンガ効果を精度よくシミュレート可能となり、ひいては、スプリングバックを精度よく解析可能となる。
【０１３１】
図７に示すように、前記折れ線は、除荷および逆負荷時における実際の応力−ひずみ関係を表すグラフにできる限り一致するように取得される。具体的には、そのグラフ上に、同一直線上に位置せず、かつ、相互に一定距離以上離れた３点を暫定的に設定することを、そのグラフのうち折れ線として抽出される可能性のある部分の一方の端点から他方の端点に向かって順次行い、それら暫定的な３点により構成される暫定的な折れ線の角度が極大となったときに、それら３点を最終的な３点に決定するとともに、それら最終的な３点により構成される折れ線を最終的な折れ線に決定する。
【０１３２】
それら３点は、除荷開始点に最も近い第１実測点６０と、次に近い第２実測点６２と、最も遠い第３実測点６４とから構成される。それら３つの実測点６０，６２，６４における各応力−ひずみ関係は、実際の応力−ひずみ関係を表すグラフ上に位置することから、実測値と一致する。
【０１３３】
なお、図７に示す例においては、折れ線を構成する２本の直線部のうち、除荷開始点（図示しない）に近い方は、前記平行四辺形の第２直線部５２と完全にではなく部分的に一致し、また、他方の直線は、平行四辺形の第４直線部５６と完全にではなく部分的に一致している。
【０１３４】
ここで、前記材料特性値取得プログラムの内容を詳述する。
【０１３５】
図８には、この材料特性値取得プログラムがフローチャートで表されている。まず、ステップＳ１０１（以下、単に「Ｓ１０１」で表す。他のステップについても同じとする）において、実験値を表すデータがＭＯ３２からコンピュータ１２のメモリ２２に取り込まれる。
【０１３６】
次に、Ｓ１０２において、予ひずみε^Ｐ _ｇに順に付与すべき番号ｉが１に設定される。その後、Ｓ１０３において、今回の予ひずみε^Ｐ _ｇ（ｉ）に関連する応力−ひずみ関係を表すデータがメモリ２２から取り込まれる。
【０１３７】
続いて、Ｓ１０４において、その取り込まれたデータに基づき、今回の予ひずみε^Ｐ _ｇ（ｉ）に関連する応力−ひずみ関係を表すグラフが想定されるとともに、そのグラフが折れ線で近似させられる。それにより、平行四辺形上の前記３つの実測点６０，６２，６４が前述のようにして抽出される。その後、Ｓ１０５において、第２実測点６２と第３実測点６４とを結ぶ第２直線５２の勾配を角度φで算出することにより、今回の予ひずみε^Ｐ _ｇ（ｉ）に対応する塑性硬化係数Ｈ（ｉ）（＝ｔａｎφ）が算出される。
【０１３８】
続いて、Ｓ１０６において、図９において式（１）で示す前述の式を用いることにより、今回の予ひずみε^Ｐ _ｇ（ｉ）に対応する材料軟化係数β（ｉ）が算出される。本実施形態においては、材料軟化係数βが０＜β＜１の範囲内となること、すなわち、材料硬化モデルが複合硬化モデルとなることを前提に算出される。以下、その算出方法を具体的に説明する。
【０１３９】
上記式（１）を用いて材料軟化係数βを算出する際、図９において式（２）で示す第１条件と、すべての予ひずみε^Ｐ _ｇ　について、式（３）で示す条件が成立するという第２条件とが制約条件として採用される。また、上記式（１）における「Ｙ_１」と「Ｈ」は、実験値により取得されたものを採用することにする。
【０１４０】
なお、再降伏応力Ｙ_１については、種々の実験により、予ひずみε^Ｐ _ｇに関して単調増加を示すことが確認されており、予ひずみε^Ｐ _ｇ＝０のときには、
Ｙ_１＝−Ｙ_０
という関係が成立し、また、すべての予ひずみε^Ｐ _ｇについて、
Ｙ_１＞−Ｙ_０
という関係が成立する。
【０１４１】
上記式（２）の第１条件を考慮し、かつ、上記式（１）を用いると、式（４）が得られる。この式（４）は式（５）に変形できる。この式（５）を用い、かつ、上記式（３）で示す第２条件を考慮すると、式（６）が得られ、この式（６）から式（７）が得られる。
【０１４２】
また、上記式（３）で示す第２条件を考慮し、かつ、上記式（１）を「β」について解くと、式（８）が得られる。
【０１４３】
そして、本実施形態においては、上記式（７）で示す条件が成立するか否かを判定し、成立しない場合には、成立するように今回の予ひずみε^Ｐ _ｇ（ｉ）が修正される。この修正は例えば、今回の予ひずみε^Ｐ _ｇ（ｉ）を、
１．０１×（−Ｙ_１／Ｈ）
として計算したり、実験値のうち上記（７）式で示す条件を満たすもののうち最も小さいものを選択することによって行われる。
【０１４４】
上記式（７）で示す条件が成立したならば、必要に応じて修正された今回の予ひずみε^Ｐ _ｇ（ｉ）と、実験値により取得された初期降伏応力Ｙ_０，再降伏応力Ｙ_１および塑性硬化係数Ｈ（ｉ）とを上記式（１）に代入することにより、今回の材料軟化係数β（ｉ）を算出する。
【０１４５】
さらに、その算出された材料軟化係数β（ｉ）が、前記式（２）で示す第１条件が成立するか否か、すなわち、０＜β＜１の範囲内にあるか否かを判定し、あればその算出された材料軟化係数β（ｉ）が最終値とされる。その第１条件が成立しない場合には、その第１条件が成立するように材料軟化係数β（ｉ）が修正され、さらに、それに合わせて初期降伏応力Ｙ_０も修正される。
【０１４６】
材料軟化係数β（ｉ）および初期降伏応力Ｙ_０の修正は例えば、もとの材料軟化係数β（ｉ）が０以下である場合には、材料軟化係数β（ｉ）を０より大きい設定値、例えば、０．１に変更するとともに、材料軟化係数β（ｉ）が０．１である場合の初期降伏応力Ｙ_０を上記式（１）を用いて計算する。その計算された初期降伏応力Ｙ_０が上記式（３）で示す条件を満たす場合には、材料軟化係数βの最新値β（ｉ）および初期降伏応力Ｙ_０がそれぞれそのまま最終値として採用される。
【０１４７】
これに対して、変更された材料軟化係数β（ｉ）の下での初期降伏応力Ｙ_０が上記式（３）で示す条件を満たさない場合には、材料軟化係数β（ｉ）が例えば０．１増加させられ、同様にして、新たな初期降伏応力Ｙ_０が計算される。このような処理は、材料軟化係数βが１を超えない範囲内で繰り返され、その結果、材料軟化係数β（ｉ）が、実験値に基づく再降伏応力Ｙ_１および塑性硬化係数Ｈ（ｉ）の下に、上記式（１）ないし式（３）のすべてを満たすように修正されることになる。
【０１４８】
以上、もとの材料軟化係数β（ｉ）が０以下である場合の修正方法を説明したが、１以上である場合にも、同様な方法で修正が行われる。
【０１４９】
図１０には、Ｓ１０６の詳細が材料軟化係数算出ルーチンとしてフローチャートで表されている。
【０１５０】
まず、Ｓ２５１において、前記抽出された第２実測点６２を用いることにより、再降伏応力Ｙ_１が決定される。次に、Ｓ２５２において、その決定された再降伏応力Ｙ_１と、前記算出された塑性硬化係数Ｈ（ｉ）と、今回の予ひずみε^Ｐ _ｇ（ｉ）とが図９の式（７）で示す条件を満たすか否かが判定される。
【０１５１】
今回は、満たすと仮定すれば、判定がＹＥＳとなり、Ｓ２５３において、上記決定された再降伏応力Ｙ_１と、前記算出された塑性硬化係数Ｈ（ｉ）と、実験値により取得された初期降伏応力Ｙ_０とを前記式（１）に代入することにより、今回の材料軟化係数β（ｉ）が算出される。その後、Ｓ２５４に移行する。
【０１５２】
これに対して、今回は、前記式（７）で示す条件を満たさないと仮定すれば、Ｓ２５２の判定がＮＯとなり、Ｓ２５５において、今回の予ひずみε^Ｐ _ｇ（ｉ）が、前述のようにして、前記式（７）で示す条件を満たすように修正される。
【０１５３】
その後、Ｓ２５３において、修正された予ひずみε^Ｐ _ｇ（ｉ）の下で前記式（１）を用いることにより今回の材料軟化係数β（ｉ）が算出される。続いて、Ｓ２５４に移行する。
【０１５４】
いずれの場合にも、その後、Ｓ２５４において、その算出された材料軟化係数β（ｉ）が、図９の式（２）で示す条件を満たすか否かが判定される。今回は、満たすと仮定すれば、判定がＹＥＳとなり、上記算出された材料軟化係数β（ｉ）と、実験値に基づく初期降伏応力Ｙ_０とがそのまま最終値として採用される。以上で本ルーチンの一回の実行が終了する。
【０１５５】
これに対して、今回は、満たさないと仮定すれば、判定がＮＯとなり、Ｓ２５６において、今回の材料軟化係数β（ｉ）が前述のようにして修正される。その後、Ｓ２５７において、修正された材料軟化係数β（ｉ）の下に、前記式（１）を用いることにより、初期降伏応力Ｙ_０が算出される。
【０１５６】
続いて、Ｓ２５８において、その算出された初期降伏応力Ｙ_０が０より大きいか否かが判定される。今回は、０より大きいと仮定すれば、判定がＹＥＳとなり、その修正された材料軟化係数β（ｉ）と、その算出された初期降伏応力Ｙ_０とがそれぞれ最終値として採用される。以上で、本ルーチンの一回の実行が終了する。
【０１５７】
また、今回は、その算出された初期降伏応力Ｙ_０が０以下であると仮定すれば、Ｓ２５８の判定がＮＯとなり、Ｓ２５６に戻る。このＳ２５６においては、材料軟化係数β（ｉ）が再度、前述のようにして修正され、その後、Ｓ２５７において、再度修正された材料軟化係数β（ｉ）の下に新たな初期降伏応力Ｙ_０が算出される。
【０１５８】
続いて、Ｓ２５８において、前記の場合と同様にして、その算出された初期降伏応力Ｙ_０が０より大きいか否かが判定される。Ｓ２５６ないしＳ２５８の実行が何回か繰り返された結果、Ｓ２５８の判定がＹＥＳとなれば、最終的に修正された材料軟化係数β（ｉ）と初期降伏応力Ｙ_０とがそれぞれ最終値として採用される。以上で本ルーチンの一回の実行が終了する。
【０１５９】
以上のようにしてＳ１０６の実行が終了すれば、その後、図８のＳ１０７において、第１実測点６０と第２実測点６２とをつなぐ第２直線５２の勾配が角度γで算出され、さらに、ｔａｎγとして、今回の予ひずみε^Ｐ _ｇ（ｉ）に対応するヤングＥ（ｉ）が算出される。その算出されたヤング率Ｅ（ｉ）から、予め定められた規則に従ってせん断弾性係数Ｇ（ｉ）が誘導される。このようにして、ヤング率Ｅ（ｉ）とせん断弾性係数Ｇ（ｉ）とが取得される。
【０１６０】
続いて、Ｓ１０８において、番号ｉが最大値ｉＭＡＸ以上であるか否かが判定される。複数の実験値のすべてが材料特性値を決定するために利用されたか否かが判定されるのである。今回は、番号ｉが最大値ｉＭＡＸ以上ではないと仮定すれば、判定がＮＯとなり、Ｓ１０９において、番号ｉが１増加させられ、その後、Ｓ１０３に移行する。以下、Ｓ１０３ないしＳ１０８が、新たな予ひずみε^Ｐ _ｇ（ｉ）に関連して実行されることになる。
【０１６１】
Ｓ１０３ないしＳ１０９の実行が何回か繰り返された結果、Ｓ１０８の判定がＹＥＳとなれば、Ｓ１１０において、予ひずみε^Ｐ _ｇの変化領域が、設定された最小値と最大値との間において、一定間隔で複数の区分に分割される。
【０１６２】
このとき、各区分に対応する予ひずみε^Ｐ _ｇに対応する材料特性値が存在しない場合には、材料特性値が存在しない予ひずみε^Ｐ _ｇに近接する複数の予ひずみε^Ｐ _ｇであって材料特性値が存在するものを内挿することにより、必要な材料特性値を生成する。その結果、予ひずみε^Ｐ _ｇと材料特性値との関係が、図１１に表形式で表されるように、予ひずみε^Ｐ _ｇの複数の代表値について離散的に取得されることになる。
【０１６３】
以上で、この材料特性値取得プログラムの一回の実行が終了する。
【０１６４】
次に、前記成形シミュレーション・プログラムの内容を詳述する。
【０１６５】
図１２には、その成形シミュレーション・プログラムの内容がフローチャートで概念的に表されており、そのうちのステップＳ４０３の詳細が図１３および図１４に要素応力計算ルーチンとしてフローチャートで表されている。以下、この成形シミュレーション・プログラムの内容を説明するが、それに先立ち、その前提となる理論および原理を説明する。
【０１６６】
この成形シミュレーション・プログラムにおいては、板材のプレス成形による変形が時間増分形の有限要素法を用いて解析される。
【０１６７】
ここに、有限要素法は、よく知られているように、材料を複数の有限要素に分割し、各要素の各節点における応力と変位とを互いに関連付ける剛性マトリクスを用い、静的な問題に対しては各要素の各節点における外力と内力とのつり合いを条件に、連立方程式である剛性方程式の繰返し計算をそれの解が真の解に収束するまで行う数値解析法である。本実施形態においては、複数の要素剛性マトリクスが集合して成る全体剛性マトリクスの各要素が未知数であって、求めるべき解とされ、その解が真の解に収束するまで行う繰返し計算が行われる。
【０１６８】
本実施形態においては、その有限要素法が増分計算法に組み合わされて実行される。増分計算法は、よく知られているように、非線形挙動や動的過程中の材料挙動を解析するために、時間を微小区間に分割し、それぞれの時間増分で計算された各物理量（変位、回転角、応力、ひずみ）の増分を順次累積して行くことによって、所定の時間までの解析を進める数値解析法である。
【０１６９】
本実施形態においては、プレス成形中のポンチの一回の動作に伴う材料挙動が、各々微小時間を有する複数の計算ステップに分割されており、それら複数の計算ステップの終了により、一連の成形シミュレーションが終了する。さらに、本実施形態においては、増分計算法として動的な陽解法が採用されている。
【０１７０】
この動的な陽解法により、ばねの運動方程式と等価な動的な運動方程式が解かれる。具体的には、各シェル要素の各節点の加速度が計算され、その加速度の時間積分によって各節点の速度が計算され、その速度の時間積分によって各節点の変位が計算される。
【０１７１】
この成形シミュレーション・プログラムにおいては、さらに、弾塑性体のための応力更新アルゴリズムに基づき、応力およびひずみが解析される。この応力更新アルゴリズムは、フォン・ミーゼス（Ｖｏｎ　Ｍｉｓｅｓ）の降伏条件のもと、各シェル要素に時刻ｔから時刻ｔ＋Δｔまでの各回の計算ステップ中に生じたひずみ増分に対する応力増分が解析される。この応力更新アルゴリズムは、図１３および図１４を参照して後述される要素応力計算ルーチンとして具体化されている。
【０１７２】
この応力更新アルゴリズムにおいては、今回の計算ステップにおける変位増分ｕと回転角増分θとが計算される。変位増分ｕは、本来、「ｄｕ」または「Δｕ」として表記すべきであるが、ここでは単に「ｕ」として表記する。また、変位増分ｕは、各成分の変位増分ｕ，ｖ，ｗの総称である。同様にして、回転角増分θは、本来、「ｄθ」または「Δθ」として表記すべきであるが、ここでは単に「θ」として表記する。また、回転角増分θは、各成分の回転角増分θｙ，θｙ，θｚの総称である。
【０１７３】
具体的には、変位増分ｕの計算は、力のつり合いを考慮した増分型の全体剛性方程式を用いることによって行われる。これに対し、回転角増分θの計算は、モーメントのつり合いを考慮した増分型の全体剛性方程式を用いることによって行われる。いずれについても、全体剛性方程式は、未知の成分を有する全体剛性マトリクスを含み、かつ、その全体剛性マトリクスの未知の成分は、今回の計算ステップにおいては前回の計算ステップまでの計算値に等しくされる。本実施形態においては、全体剛性マトリクスの未知の成分が有限要素法によって求められる。
【０１７４】
前記ひずみ増分は、等方ひずみ増分ε_ｍと偏差ひずみ増分ｅ_ｉｊとに分離される。その等方ひずみ増分ε_ｍは、各成分ごとのひずみ増分ε_ｉｊを用いて計算され、そのひずみ増分ε_ｉｊは、後述のひずみ近似関数を用いて計算される。
【０１７５】
等方ひずみ増分ε_ｍから静水圧増分３Ｋε_ｍ（Ｋ：体積弾性係数）が求められる。この静水圧増分３Ｋε_ｍに前回の計算ステップまでの静水圧成分σ_ｉｊ ^０を加えることにより、静水圧成分σ_ｍが更新される。すなわち、更新された静水圧成分σ_ｍ ^ｔ＋Δｔが計算されるのであり、このことが図１５に式（１１）で示されている。
【０１７６】
各成分ごとのひずみ増分ε_ｉｊから等方ひずみ増分ε_ｍを差し引くことにより、各成分ごとの偏差ひずみ増分ｅ_ｉｊが求められる。
【０１７７】
この偏差ひずみ増分ｅ_ｉｊから偏差応力成分２Ｇｅ_ｉｊ（Ｇ：せん断弾性係数）が求められる。したがって、偏差応力成分２Ｇｅ_ｉｊは、平面応力状態において弾性域内において求められることになる。この偏差応力成分２Ｇｅ_ｉｊに前回の計算ステップまでの偏差応力成分Ｓ_ｉｊ ^０を加えることにより、試行応力（偏差応力Ｓの試行値）Ｓ_ｉｊ ^{Ｔｒｉａｌ}が求められる。すなわち、更新された試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}が計算されるのであり、このことが式（１２）で示されている。
【０１７８】
この試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}を式（１３）に代入することにより、相当応力の試行解σ^{Ｔｒｉａｌ}が求められる。
【０１７９】
この相当応力の試行解σ^{Ｔｒｉａｌ}が弾性域内にあり、各シェル要素の材料が未だ降伏していない場合には、上記試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}がそのまま、更新された偏差応力成分（今回の偏差応力成分）Ｓ_ｉｊ ^{Ｆｉｎａｌ}とされる。
【０１８０】
一方、応力空間を表す図１６の（ａ）に示すように、相当応力の試行解σ^{Ｔｒｉａｌ}が塑性域に達すると、各シェル要素の材料が降伏したと判断される。この場合には、同図の（ｂ）に示すように、図１５に式（１４）で示すラジアルリターンアルゴリズム（半径引戻し解法）により、試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}が降伏曲面上へ引き戻される。
【０１８１】
そのラジアルリターンアルゴリズムにおいては、相当塑性ひずみ増分ε^Ｐから、降伏曲面の半径を表す相当応力σと、最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とが計算される。さらに、それら両計算値が比較され、相当応力σに一致するときの最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値が最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の最終値に決定される。
【０１８２】
なお付言すれば、前述のフォン・ミーゼスの降伏条件を採用する場合には、例えば、図１５の式（１６）を用いることより増分Δε^Ｐを計算し、その後、式（１４）を用いることにより最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}を計算することが可能である。
【０１８３】
各シェル要素の材料が降伏した場合には、その引き戻された最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の最終値が、更新された偏差応力成分として求められる。
【０１８４】
そして、図１５に式（１５）で示すように、上記更新された静水圧成分σ_ｍ ^ｔ＋Δｔと、上記更新された偏差応力成分Ｓ_ｉｊ ^{Ｆｉｎａｌ}とにより、各成分ごとに更新された応力σ_ｉｊ ^ｔ＋Δｔが計算される。
【０１８５】
以上説明した応力更新の１回の計算ステップは、板材を構成するすべてのシェル要素について順に実行される。
【０１８６】
従来のシェル要素では、曲げ変形の前後で断面の平面性が維持されると仮定されており、そのため、シェル要素の中立面は常に板厚中央面と一致するとともに、断面回転角θｙが、曲げ方向位置を表すｘ座標値には依存するが、板厚方向位置を表すｚ座標値には依存しない値とされる。
【０１８７】
したがって、断面回転角θｙは、図１７の式（２１）で表され、曲げ方向の垂直ひずみεｘは式（２２）で表され、板厚方向のせん断ひずみγｘｚは式（２３）で表される。ただし、それらの式において、「ｕ」は、ｘ方向の変位を表し、「ｗ」は、ｚ方向の変位を表している。
【０１８８】
よって、従来のシェル要素を用いる場合には、図１８に示すように、垂直ひずみεｘは板厚方向に線形に分布する一方、せん断ひずみγｘｚが板厚方向に沿って変化しないように分布する。
【０１８９】
しかし、本発明者らの研究により、そのような仮定は、曲げ変形させられるべき板材が薄板である場合には妥当であるが、厚板である場合のように、板材の曲げの程度が厳しい場合には妥当性を失うことが判明した。
【０１９０】
図１９には、断面回転角θｙと、板材内の、それの曲げ内面からの距離ｄｓとの関係が曲げの厳しさＬｂによって変化する様子がグラフで表されている。このグラフは、本発明者らの実験結果に基づいて作成されたものである。ここに、曲げの厳しさＬｂは、曲げ半径Ｒを板厚ｔで割り算した値である曲げ指標Ｒ／ｔを用い、その値が小さいほど曲げの厳しさが強くなるように表現される。
【０１９１】
図１９のグラフによれば、曲げの厳しさＬｂが３以上である領域（曲げの厳しさが弱い領域であって、同図においては「Ｌｂ＝３」という表記に関連付けられた水平線の下側の領域）においては、板厚方向位置にかかわらず断面回転角θｙがほとんど変化せず、よって、従来のシェル要素についての仮定が妥当であることが分かる。しかし、曲げの厳しさＬｂが３より小さい領域（曲げの厳しさが強い領域であって、上記水平線の上側の領域）においては、板厚方向位置によって断面回転角θｙが顕著に変化し、よって、従来のシェル要素についての仮定が妥当ではないことが分かる。
【０１９２】
そこで、本発明者らは、厚板中の任意の部位の断面回転角θｙをシミュレーションによって計算した。このシミュレーションにおいては、厚板が複数のソリッド要素に仮想的に分割される。ソリッド要素モデルを用いてシミュレーションが行われるのである。
【０１９３】
このシミュレーションのために実行されるプログラムが前記成形シミュレーション・プログラムに組み込まれている。図２０には、そのプログラムがひずみ近似関数特定ルーチンとして概念的に表されている。
【０１９４】
このひずみ近似関数特定ルーチンにおいては、まず、Ｓ６０１において、曲げ指標Ｒ／ｔが異なる複数の代表板材（潜在的な等価板材）のそれぞれにつき、有限要素法により、各ソリッド要素ごとに変位ｕ（各成分の変位ｕ，ｖ，ｗ）と回転角θ（各成分の回転角θｙ，θｙ，θｚ）とが計算される。次に、Ｓ６０２において、それら計算された変位ｕと回転角θとに基づき、かつ、図２５の式（３２）および式（３３）を用いることにより、ひずみがｚ座標値に関連付けて計算される。
【０１９５】
続いて、Ｓ６０３において、その計算されたひずみとｚ座標値との関係に基づき、ひずみの板厚方向における分布を近似的に定義するひずみ近似関数が特定される。その後、Ｓ６０４において、その特定されたひずみ近似関数が、各板材の曲げ指標Ｒ／ｔに関連付けられる。以上で、このひずみ近似関数特定ルーチンの一回の実行が終了する。
【０１９６】
このひずみ近似関数特定ルーチンは、スプリングバックを解析すべき板材ごとに実行することは不可欠ではない。このひずみ近似関数特定ルーチンの過去の実行により、実質的に同じ曲げの厳しさＬｂに関してひずみ近似関数が既に特定されている場合には、その特定されたひずみ近似関数を利用することにより、ひずみ近似関数特定ルーチンの新たな実行を省略することができるからである。
【０１９７】
本発明者らは、ソリッド要素モデルを用いたシミュレーションの精度を検証するため、厚板につき、実験によって断面回転角θｙと曲げ内面からの距離ｄｓとの関係を測定した。さらに、本発明者らは、その測定された関係を表す実験値と上述のシミュレーションによる計算値とを互いに比較した。
【０１９８】
図２１には、厚板につき、断面回転角θｙと距離ｄｓとの関係が実験値と計算値との双方によって示されている。同図から明らかなように、それら実験値と計算値との間には良好な一致が認められる。したがって、厚板につき、実験値に代えて計算値を用いることが可能であり、よって、実験なしでも、断面回転角θｙと距離ｄｓとの関係、すなわち、断面回転角θｙのｚ座標値への依存性を取得することが可能である。
【０１９９】
図２２には、上記シミュレーションによる計算値から計算された板厚方向のせん断ひずみγｘｚ（前記ひずみの一例）と曲げ内面からの距離ｄｓとの関係が曲げ指標Ｒ／ｔの各値に関連付けてグラフで示されている。同図には、さらに、曲げ指標Ｒ／ｔの値がそれぞれ１．８３，１．６１，２．０５である場合のせん断ひずみγｘｚと距離ｄｓとの関係を近似するひずみ近似関数がグラフで示されている。このひずみ近似関数は、距離ｄｓを表す変数ｚを変数とする関数によってせん断ひずみγｘｚを定義する式である。このひずみ近似関数は、例えば、２次関数、３次関数、４次関数または５次関数である多項式として構成することが可能である。
【０２００】
このひずみ近似関数を利用すれば、シェル要素を用いた成形シミュレーションであっても、ソリッド要素を用いた場合とほぼ同等の精度で、ひずみの板厚方向における分布を計算することが可能となる。そして、図２３には、せん断ひずみγｘｚに関し、ソリッド要素を用いた成形シミュレーションによる計算値と、本実施形態のシェル要素を用いた成形シミュレーションによる計算値とがグラフで示されている。両者は良好な一致を示している。
【０２０１】
図２３には、さらに、比較のために、従来例のシェル要素を用いた成形シミュレーションによる計算値も示されている。従来例のシェル要素においては、せん断ひずみγｘｚがｚ座標値に依存しない一定値であると仮定されている。
【０２０２】
図２４には、本実施形態のシェル要素を用いて計算されたせん断ひずみγｘｚの分布と、従来例のシェル要素を用いて計算されたせん断ひずみγｘｚの分布とがグラフで表されている。本実施形態のシェル要素を用いて計算されたせん断ひずみγｘｚの分布は、前記ひずみ近似関数が２次関数である場合の一例を示している。このように、本実施形態のシェル要素を用いる場合には、従来例のシェル要素を用いる場合に比較し、厚板の材料特性の実際、すなわち、せん断ひずみγｘｚの分布の実際が精度よく再現される。
【０２０３】
従来例のシェル要素とは異なり、本実施形態のシェル要素については、断面回転角θｙがｘ座標値のみならずｚ座標値にも依存するが、このことが図２５において式（３１）で表されている。
【０２０４】
このような依存性により、本実施形態のシェル要素については、曲げ方向の垂直ひずみεｘが図２５の式（３２）で表され、板厚方向のせん断ひずみγｘｚは式（３３）で表される。
【０２０５】
図２６には、本実施形態のシェル要素を用いた計算値が、せん断ひずみγｘｚと垂直ひずみεｘとのそれぞれの、板厚方向における分布に関して示されている。せん断ひずみγｘｚの分布は、前記ひずみ近似関数を２次より高次の関数である場合の一例を示している。本実施形態のシェル要素によれば断面回転角θｙをｚ座標値に依存させることが可能になり、その結果、せん断ひずみγｘｚのみならず垂直ひずみεｘについても実際の分布が精度よく再現される。
【０２０６】
以上、本実施形態のシェル要素と従来例のシェル要素とをひずみに関して互いに比較したが、次に、板厚方向の垂直応力σｚに関して互いに比較する。
【０２０７】
前述のように、従来例のシェル要素においては、図２７に示すように、垂直応力σｚが無視されて０であると仮定される。これに対し、本実施形態のシェル要素においては、垂直応力σｚが考慮される。
【０２０８】
板厚方向の垂直応力σｚは、厚板の曲げ変形による寄与分σｚｓと、厚板がポンチから受ける接触圧ｐによる寄与分σｚｃとの和であると考えることができる。寄与分σｚｃは、接触圧ｐによって厚板に生ずるせん断変形に起因する。
【０２０９】
垂直応力σｚは、厚板（厚肉円筒）内の微小要素のｚ方向（厚肉円筒の半径方向）における力のつり合いに着目することにより、垂直応力σｘから解析することが可能である。その力のつり合いは、極座標で記述するとともに時間的に離散化すれば、図２８における式（４１）で表される。この式（４１）における「ｒ」は、ｚ座標軸上の変数であり、「ｚ」に置き換えることが可能である。
【０２１０】
本実施形態においては、まず、垂直応力σｚを無視して垂直応力σｘが暫定的に解析され、その暫定値に基づき、かつ、上記式（４１）を用いて、垂直応力σｚの板厚方向における分布を近似的に表す応力近似関数が特定される。垂直応力σｘの暫定値は、厚板の平面応力状態における垂直応力σｘに一致する。
【０２１１】
その応力近似関数の特定は具体的に次のようにして行うことが可能である。
【０２１２】
まず、応力近似関数の基本形を選択する。この基本形は、ｎ次関数により構成したり、指数関数により構成することが可能である。この基本形は、厚板の曲げ内面がポンチから受ける接触圧を考慮するか否かによって異なる。
【０２１３】
図２８の式（４２）は、ポンチからの接触圧を考慮しない応力近似関数の一例である。この式（４２）で表される応力近似関数において「ａ」は、厚板の曲げ内面のｚ方向位置を表し、また、「ｂ」は、厚板の曲げ外面のｚ方向位置を表している。中立面のｚ方向位置をｒ_０、板厚をｔ_０でそれぞれ表すと、「ａ」は式（４３）で、「ｂ」は式（４４）でそれぞれ表される。
【０２１４】
したがって、上記式（４２）で表される垂直応力σｚは、厚板の両表面上において０となり、このことは、厚板の曲げ内面がポンチから受ける接触圧が無視されることを示している。その接触圧を考慮して応力近似関数を特定する手法については後述する。
【０２１５】
垂直応力σｚを式（４２）で表される応力近似関数ｆ（ｚ）で表すこととすれば、前記式（４１）が式（４５）に変形される。
【０２１６】
上記式（４２）は式（４６）で表すように展開できる。また、応力近似関数ｆ（ｚ）の一次導関数ｆ’（ｚ）は、式（４７）で表される。それら式（４６）および（４７）を上記式（４５）に代入すれば、式（４８）が得られる。
【０２１７】
式（４８）中の各係数に式（４３）および式（４４）を代入するとともに、変数ａ０，ａ１，ｌ，ｍを用いた変数置換を行うと、式（４８）中の各係数が図２８の式（４９）で示すように整理される。この整理の結果を用いれば、図２８の式（４８）が図２９の式（５０）に変形される。
【０２１８】
式（５１）で示す変数置換を行うと、式（５０）が式（５２）に変形される。この式（５２）においては、「Ｘ」と「Ｙ」と「Ｗ」と「Ｚ」とが既知である一方、「ｌ」と「ｍ」と「Ａ’」とが未知である。
【０２１９】
ところで、本実施形態においては、各シェル要素が、それの板厚方向に積層された複数のレイヤに関連付けられるとともに、各シェル要素ごとに、垂直応力σｘの板厚方向における分布を表す関数（各レイヤの位置すなわちｚ座標値を変数とする）が予め定義される。それにより、本実施形態においては、各シェル要素につき、各ｚ座標値ごとに垂直応力σｘが求められ、それらｚ座標値と垂直応力σｘとの関係を前提にして最小２乗法が式（５２）に適用されることにより、未知数ｌ，ｍ，Ａ’が計算される。
【０２２０】
図３０の式（５３）は、図２９の式（５２）について２乗和を表している。この式（５３）を未知数ｌ，ｍ，Ａ’に関してそれぞれ偏微分すれば、式（５４）が得られ、これは式（５５）に変形できる。この式（５５）で表される３つの連立方程式を解けば、未知数ｌ，ｍ，Ａ’が求められ、これにより、応力近似関数ｆ（ｚ）が特定される。
【０２２１】
なお付言すれば、関数によって板厚方向分布が定義される変数に、垂直応力σｚを選んだり、せん断応力τｘｚを選んだり、両者を選ぶことが可能である。
【０２２２】
以上、応力近似関数ｆ（ｚ）の例として４次関数である場合を説明したが、次に、３次関数である場合を説明する。
【０２２３】
垂直応力σｚを図３１の式（５６）で表される応力近似関数ｆ（ｚ）で表すこととすれば、応力近似関数ｆ（ｚ）の一次導関数ｆ’（ｚ）は、式（５７）で表される。それら式（５６）および（５７）を前記式（４５）に代入すれば、式（５８）が得られる。この式（５８）を変数ｃに関して整理すれば、式（５９）が得られる。
【０２２４】
式（６０）で示す変数置換を行うと、上記式（５９）が式（６１）に変形される。この式（６１）においては、「Ｘ」と「Ｙ」と「Ｚ」とが既知である一方、「ｃ」と「Ａ’」とが未知である。各シェル要素につき、各レイヤのｚ座標値と各レイヤの垂直応力σｘとの関係を前提にして最小２乗法が式（６１）に適用されることにより、未知数ｃ，Ａ’が計算される。式（６１）の２乗和が図３２の式（６２）で表されている。
【０２２５】
式（６２）を未知数ｃに関して偏微分すれば、式（６３）が得られ、これは式（６４）に変形できる。これに対し、式（６２）を未知数Ａ’に関して偏微分すれば、式（６５）が得られ、これは式（６６）に変形できる。それら式（６４）および（６６）の連立方程式を解くことにより、未知数ｃ，Ａ’が求められ、これにより、応力近似関数ｆ（ｚ）が特定される。
【０２２６】
以上、垂直応力σｚの近似関数をポンチからの接触圧を考慮しないで特定する手法を説明したが、その接触圧を考慮する場合には、その近似関数の基本形が例えば、図３３において式（６７）で表される。その場合には、厚板の曲げ内面（ｚ＝ａ）で垂直応力σｚが接触圧ｐに等しく、一方、曲げ外面（ｚ＝ｂ）で垂直応力σｚが０となる。このことが式（６８）で表されている。
【０２２７】
接触圧ｐを考慮して応力近似関数ｆ（ｚ）を求める別の例を説明する。この例においては、垂直応力σｚが、図３４の式（６９）に示すように、前述の、接触圧ｐの寄与分σｚｃと曲げ変形の寄与分σｚｓとの和として表現される。寄与分σｚｃは、例えば、式（７０）に示す３次関数で定式化される。この寄与分σｚｃは、曲げ内面においては接触圧ｐに等しく、曲げ外面においては０であり、このことが式（７１）で表されている。これに対し、寄与分σｚｓは、式（７２）で表現される。
【０２２８】
式（７０）を変数ｒ（変数ｚと等しい）に関して偏微分すると、式（７３）が得られる。したがって、この場合、微小円筒要素における板厚方向（半径方向）の力のつり合いは式（７４）で表される。この式（７４）を前述と同様な手法で解くことにより、垂直応力σｚを近似する応力近似関数ｆ（ｚ）が求められる。
【０２２９】
垂直応力σｚを求めるためにポンチからの接触圧ｐを考慮しない場合であっても考慮する場合であっても、その後、前述のように、その垂直応力σｚに基づいて垂直応力σｘの暫定値（平面応力状態での値）が修正される。以下、この修正手法を説明する。
【０２３０】
プラントル−ルイス（Ｐｒａｎｄｔｌ−Ｒｅｕｓｓ）の流れ法則によれば、各成分の偏差応力Ｓについて図３５の式（８１）が成立する。また、垂直応力σｚを無視する場合には、各成分の偏差応力ｓについて式（８２）が成立する。これに対し、垂直応力σｚを考慮する場合には、各成分の偏差応力ｓについて式（８３）が成立する。式（８２）における「σｘ」および「σｙ」は、垂直応力σｚを考慮しない場合のｘ方向垂直応力およびｙ方向垂直応力をそれぞれ意味し、一方、式（８３）における「σｘ’」および「σｙ’」は、垂直応力σｚを考慮する場合のｘ方向垂直応力およびｙ方向垂直応力をそれぞれ意味している。
【０２３１】
式（８２）中の第１式の右辺の値と式（８３）中の第１式の右辺の値とが互いに一致することから、式（８４）が誘導される。また、式（８２）中の第２式の右辺の値と式（８３）中の第２式の右辺の値とが互いに一致することから、式（８５）が誘導される。
【０２３２】
それら式（８４）および式（８５）から、垂直応力σｘ’については式（８６）、垂直応力σｙ’については式（８７）がそれぞれ誘導される。それら式（８６）および式（８７）によれば、垂直応力σｚを考慮する垂直応力σｘ’，σｙ’は、考慮しない垂直応力σｘ，σｙに対して垂直応力σｚの分だけ増加することが分かる。この増加は、厚板の横断面全体に作用する垂直応力σｘ’，σｙ’の各総和である断面力（垂直応力σｚを考慮した断面力）が、垂直応力σｘ，σｙの各総和である断面力（平面応力状態での断面力）に対して垂直応力σｚの総和の分だけ増加することを意味する。
【０２３３】
そこで、本実施形態においては、垂直応力σｚを考慮した断面力が、平面応力状態での断面力に維持されるように、中立面の板厚方向における位置ｒ_０が変更される。具体的には、垂直応力σｚを考慮した断面力が、垂直応力σの増加分だけ減少するように、中立面の板厚中央位置からの移動量が決定される。
【０２３４】
図３６には、中立面の移動量を決定する手法の一例が、その決定を垂直応力σｘのみに関して行う場合を例にとり、グラフで示されている。同図には、移動前の中立面に対応する垂直応力σｘが実線グラフ、移動後の中立面に対応する垂直応力σｘが二点鎖線グラフでそれぞれ示されている。
【０２３５】
図３６において実線グラフを、中立面が曲げ内面に接近する向きに距離Ｌ１だけ移動させると、引張側（＋側）においては、垂直応力σｚの曲げ外面における値と距離Ｌ１との積に等しい分だけ、垂直応力σｘの総和である断面力が増加する一方、圧縮側（−側）においては、垂直応力σｚの曲げ内面における値と距離Ｌ１との積に等しい分、垂直応力σｘの総和である断面力が減少する。引張側の増加分は、圧縮側の減少分に相当すると考えることができる。
【０２３６】
したがって、板材の最外表面上の応力の関数を求め、かつ、予め計算された垂直応力σｘの増加量を利用すれば、距離Ｌ１が一義的に決まり、これにより、最終的な中立面の位置が決定される。
【０２３７】
図３６に示す例においては、垂直応力σｚの符号が負である場合の一例であり、そのため、中立面の位置が、垂直応力σｘの総和である断面力が増加する向きに移動させられる。すなわち、この例においては、垂直応力σｚを考慮することによる垂直応力σｘの総和である断面力の減少分が中立面の移動によって補われ、その結果、垂直応力σｚを考慮する場合と考慮しない場合とで、垂直応力σｘの総和である断面力が変化しないようになっている。これにより、同じ節点を共有する２つのシェル要素間において断面力が、垂直応力σｚを考慮したことが原因で不連続となってしまう不自然な事態が回避される。
【０２３８】
図３７には、曲げ方向応力σｘについての計算値がグラフで表されている。それら計算値は、
（１）シェル要素モデルを用いた計算値であって、板厚方向応力σｚを考慮しない曲げ方向応力σｘを表すものと、
（２）シェル要素モデルを用いた計算値であって、板厚方向応力σｚを考慮する曲げ方向応力σｘと、
（３）シェル要素モデルを用いた計算値であって、板厚方向応力σｚを考慮するとともに力のつり合いを考慮する曲げ方向応力σｘと
を含んでいる。
【０２３９】
図３８には、さらに、ソリッド要素モデルを用いて解析された板厚方向応力σｚの値もグラフで表されている。
【０２４０】
図３８には、曲げ方向応力σｘについてのいくつかの計算値がグラフで表されている。それら計算値は、本実施形態のシェル要素を用いた計算値と、従来例のシェル要素を用いた計算値とを含んでいる。同図には、さらに、比較のため、ソリッド要素を用いて計算値もグラフで表されている。この計算値は、曲げ方向応力σの実測値と等価であると考えることができる。したがって、同図によれば、本実施形態のシェル要素を用いて計算された曲げ方向応力σｘが、従来例のシェル要素を用いて計算された曲げ方向応力σｘより精度よく実測値を反映することが分かる。
【０２４１】
以上、図１２の成形シミュレーション・プログラムの内容を概略的に説明したが、以下、具体的に説明する。
【０２４２】
この成形シミュレーション・プログラムにおいては、まず、Ｓ４０１において、成形の解析に必要なデータがＭＯ３２から読み込まれるとともに、解析の初期化が行われる。
【０２４３】
次に、Ｓ４０２において、板材に作用する荷重に関する境界条件が適用される。本実施形態においては、パッドから板材に作用する荷重に関する境界条件が適用される。
【０２４４】
続いて、Ｓ４０３において、板材を構成する複数のシェル要素のうち今回の実行対象である注目シェル要素に作用する応力が計算される。このＳ４０３の詳細が図１３および図１４に要素応力計算ルーチンとしてフローチャートで表されている。この要素応力計算ルーチンは、前記応力更新アルゴリズムを具体化するために設計されている。
【０２４５】
この要素応力計算ルーチンにおいては、まず、Ｓ５０１において、前述のように、増分型の剛性方程式において全体剛性マトリクスをそれの各成分が前回値に等しい状態で用いることにより、各シェル要素の各節点における変位増分ｕおよび回転角増分θが計算される。これらの値をもとに、板厚中央面でのひずみが計算される。次に、Ｓ５０２において、前記ひずみ近似関数に基づき、各成分ごとにひずみ増分ε_ｉｊが計算される。ここに、ひずみ増分ε_ｉｊは、垂直ひずみとせん断ひずみとの双方を含んでいる。
【０２４６】
その後、Ｓ５０３において、その計算されたひずみ増分ε_ｉｊに基づき、等方ひずみ増分ε_ｍが計算される。続いて、Ｓ５０４において、その計算された等方ひずみ増分ε_ｍに基づき、静水圧増分３Ｋε_ｍが計算される。その後、Ｓ５０５において、前回の計算ステップにおける各成分ごとの応力σ_ｉｊ ^０と、静水圧増分３Ｋε_ｍとに基づき、前述のようにして、更新された静水圧σ_ｍ ^ｔ＋ΔＴが計算される。
【０２４７】
続いて、Ｓ５０６において、前記計算された各成分ごとのひずみ増分ε_ｉｊと等方ひずみ増分ε_ｍとに基づき、各成分ごとに偏差ひずみ増分ｅ_ｉｊが計算される。続いて、Ｓ５０７において、平面応力状態で、かつ、弾性域内において、偏差応力増分２Ｇｅ_ｉｊが計算される。
【０２４８】
その後、Ｓ５０８において、前回の計算ステップにおける各成分ごとの偏差応力Ｓ_ｉｊ ^０と、偏差応力増分２Ｇｅ_ｉｊとに基づき、前述のようにして、更新された試行応力Ｓ_ｉｊ ^ｔ＋ΔＴが計算される。続いて、Ｓ５０９において、試行応力Ｓ_ｉｊ ^ｔ＋ΔＴの計算値に基づき、前述のようにして、相当応力の試行解σ^{Ｔｒｉａｌ}が計算される。
【０２４９】
その後、図１４のＳ５１０において、板材がそれの両面の少なくとも一方において降伏したか否かが判定される。今回は、降伏したと仮定すれば、判定がＹＥＳとなり、Ｓ５１１以下のステップに移行する。
【０２５０】
Ｓ５１１においては、図１５の式（１６）を用いることにより、相当塑性ひずみ増分ε^Ｐの増分Δε^Ｐが計算される。その後、Ｓ５１２において、その計算された増分Δε^Ｐに対応する相当応力σが対応相当応力σとされるとともに、その対応相当応力σと試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}とに対応する最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}が、図１５の式（１４）を用いて計算される。本実施形態においては、その最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}が前記（６）項における「暫定値」の一例である。
【０２５１】
なお付言すれば、前記フォン・ミーゼスの降伏条件を採用しない場合には、前記ラジアルリターンアルゴリズムを、Ｓ５１１およびＳ５１２の反復的実行により実現することが可能である。
【０２５２】
例えば、Ｓ５１１の各回の実行時には、相当塑性ひずみ増分ε^Ｐが一定の微小増分Δε^Ｐで増加させられる。これに対し、Ｓ５１２の各回の実行時には、相当塑性ひずみ増分ε^Ｐの暫定値（Ｓ５１０の実行によって計算される）に対応する相当応力σが対応相当応力σとして計算される。さらに、その対応相当応力σと試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}とに対応する最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値が計算される。
【０２５３】
このＳ５１２においては、さらにまた、対応相当応力σと最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とが互いに十分に一致するか否かが判定される。十分に一致する場合には、相当塑性ひずみ増分ε^Ｐの暫定値と最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とがそれぞれ最終値とされる。最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の最終値は、更新された偏差応力成分を意味する。これに対し、対応相当応力σと最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とが互いに十分に一致しない場合には、Ｓ５１１に戻り、相当塑性ひずみ増分ε^Ｐの暫定値が一定の微小増分Δε^Ｐで増加させられる。
【０２５４】
いずれの場合にも、前記ラジアルリターンアルゴリズムによって相当塑性ひずみ増分ε^Ｐの最終値と、更新された偏差応力成分Ｓ_ｉｊ ^{Ｆｉｎａｌ}とが取得されたならば、その後、Ｓ５１３において、更新された偏差応力成分Ｓ_ｉｊ ^{Ｆｉｎａｌ}と、前記更新された静水圧成分σ_ｍ ^ｔ＋Δｔとに基づき、図１５の式（１５）を用いることにより、各成分ごとに更新された応力σ_ｉｊ ^ｔ＋Δｔが計算される。
【０２５５】
以上、板材がそれの両面の少なくとも一方において降伏した場合を説明したが、いずれの面においても降伏しない場合には、Ｓ５１０の判定がＮＯとなり、Ｓ５１４に移行する。このＳ５１４においては、前回の計算ステップまでに更新された試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}と、前記更新された静水圧成分σ_ｍ ^ｔ＋Δｔとに基づき、図１５の式（１７）を用いることにより、各成分ごとに更新された応力σ_ｉｊ ^ｔ＋Δｔが計算される。
【０２５６】
いずれの場合にも、その後、Ｓ５１５において、今回の注目シェル要素の横断面に作用する曲げ方向応力σｘ（平面応力状態での値）の総和と、断面方向応力σｙ（平面応力状態での値）の総和とがそれぞれ断面力として計算される。それら断面力は、板厚方向応力σｚが０であるとして計算された値を有する。
【０２５７】
続いて、Ｓ５１６において、前述の最小２乗法により、板厚方向応力σｚの応力近似関数が特定される。その後、Ｓ５１７において、その特定された応力近似関数を満たす曲げ方向応力σｘと断面方向応力σｙとが計算される。本実施形態においては、それら計算された曲げ方向応力σｘと断面方向応力σｙとがそれぞれ前記（６）項における「仮面内応力」の一例である。
【０２５８】
なお付言すれば、この要素応力計算ルーチンは、シェル要素を用いた解析を行うにもかかわらず、板厚方向応力σｚを考慮して解析を行う機能を予め付加されている。そして、本実施形態においては、まず、板厚方向応力σｚを０と仮定して曲げ方向応力σｘと断面方向応力σｙとを求め、次に、それら曲げ方向応力σｘと断面方向応力σｙとから、前記応力近似関数を用いることにより、板厚方向応力σｚを求める。その後、上記機能を利用することにより、その求められた板厚方向応力σｚを前提にして、すなわち、板厚方向応力σｚを考慮して、曲げ方向応力σｘと断面方向応力σｙとを再度求める。
【０２５９】
続いて、Ｓ５１８において、前述のように、その計算された各応力σｘ，σｙ（垂直応力σｚを考慮した場合の値）の総和が、前記計算された断面力と一致するように、中立面が移動させられる。具体的には、中立面位置ｒ_０が補正される。この補正値は、図２５の式（３２）および式（３３）を用いたひずみの計算、応力近似関数の特定、ひずみ近似関数の特定、スプリングバックの解析等に影響を及ぼす。
【０２６０】
すなわち、本実施形態においては、その補正された中立面位置ｒ_０により規定される曲げ方向応力σｘと断面方向応力σｙとがそれぞれ前記（６）項における「最終値」の一例を構成しているのである。また、図１２におけるＳ４０９において一連の解析が終了したと判定されたときにおけるそれら曲げ方向応力σｘと断面方向応力σｙとの値がそれぞれ前記残留応力を意味している。その残留応力は、各シェル要素に関連付けてメモリ２２に記憶される。
【０２６１】
その後、Ｓ５１９において、前回の計算ステップまでの各成分のひずみ増分と、今回の偏差ひずみ増分とに基づき、図１５の式（１９）を用いることにより、更新された塑性ひずみε^Ｐ _ｉｊ ^ｔ＋Δｔが計算される。この式（１９）における「Δλ」は、式（１８）および式（１９）から求められる。図１２におけるＳ４０９において一連の解析が終了したと判定されたときにおけるその塑性ひずみε^Ｐ _ｉｊ ^ｔ＋Δｔが前記残留ひずみを意味している。その残留ひずみは、各シェル要素に関連付けてメモリ２２に記憶される。
【０２６２】
続いて、Ｓ５２０において、その更新された塑性ひずみε^Ｐ _ｉｊ ^ｔ＋Δｔに基づき、今回の注目シェル要素の板厚ｔ_０が更新される。この更新値は、応力近似関数の特定、ひずみ近似関数の特定、スプリングバックの解析等に影響を及ぼす。
【０２６３】
その後、Ｓ５２１において、Ｓ５０１ないしＳ５２０の実行が板材のすべてのシェル要素について行われたか否かが判定される。Ｓ５０１ないしＳ５２０の実行がすべてのシェル要素について行われてはいないと仮定すれば、判定がＮＯとなり、Ｓ５２２において、次のシェル要素が新たな注目シェル要素とされた後、Ｓ５０１に戻る。これに対し、Ｓ５０１ないしＳ５２０の実行がすべてのシェル要素について行われたと仮定すれば、Ｓ５２１の判定がＹＥＳとなり、Ｓ５２３において、板材を近似するシェル要素モデルに対して有限要素法が実行される。
【０２６４】
具体的には、各シェル要素ごとの要素剛性マトリクスによって全体剛性マトリクスが構築され、各節点における外力と内力とのつり合いを考慮しつつ、連立の剛性方程式が解かれる。
【０２６５】
Ｓ５２４においては、その剛性方程式の解が収束したか否かが判定される。収束しない場合には、判定がＮＯとなり、図１３のＳ５０１に戻り、新たな全体剛性マトリクスを用いてＳ５０１ないしＳ５２０が実行される。これに対し、上記剛性方程式の解が収束した場合には、Ｓ５２４の判定がＹＥＳとなり、以上で、この要素応力計算ルーチンの一回の実行が終了する。
【０２６６】
その後、図１２のＳ４０４において、前記境界条件のもと、板材の変形とポンチの形状とを考慮することにより、ポンチから板材に作用する接触力が計算される。
【０２６７】
続いて、Ｓ４０５において、ポンチの加工速度に関する境界条件を考慮することにより、各シェル要素の各節点の加速度が更新される。その後、Ｓ４０６において、ポンチと板材との接触に関する解析条件を考慮することにより、各シェル要素の各節点の加速度が更新される。続いて、Ｓ４０７において、それら節点加速度に対して時間積分を行うことにより、各節点の速度が更新され、さらに、その計算された速度に対して時間積分を行うことにより、各節点の変位が更新される。
【０２６８】
その後、Ｓ４０８において、解析結果が出力されるとともに、必要なポスト処理が行われる。続いて、Ｓ４０９において、今回の解析が終了したか否かが判定される。Ｓ４０２ないしＳ４１０から成るループの一回の実行は、増分計算法における１回の時間増分（１回の計算ステップ）に対応しており、現象的には、ポンチがダイス内に微小距離進入する動作に対応する。
【０２６９】
したがって、Ｓ４０９においては、上記ループの初回の実行開始時刻からの経過時間が、そのループの所定実行回数に対応する長さに達したか否かが判定されることにより、今回の解析が終了したか否かが判定される。現在時刻ｔにおいては、未だ解析が終了してはいないと仮定すれば、判定がＮＯとなり、Ｓ４１０において、更新された現在時刻ｔ＋Δｔ（Δｔ：時間増分）が計算される。その後、Ｓ４０２に戻る。
【０２７０】
Ｓ４０２ないしＳ４１０から成るループが何回か繰り返された結果、今回の解析が終了すれば、Ｓ４０９の判定がＹＥＳとなり、この成形シミュレーション・プログラムの一回の実行が終了する。
【０２７１】
次に、前記材料特性値決定プログラムの内容を詳述する。
【０２７２】
図３９には、この材料特性値決定プログラムがフローチャートで示されている。まず、Ｓ２０１において、材料が仮想的に分割された複数のシェル要素に順に付与された番号ｊが１に設定される。次に、Ｓ２０２において、メモリ２２から、今回のシェル要素ｊに対応する初期ひずみε^Ｐ _０（ｊ）が読み出される。
【０２７３】
その後、Ｓ２０３において、その今回の初期ひずみε^Ｐ _０（ｊ）が、前記複数の予ひずみ区分のいずれに属するかが判定され、さらに、その属すると判定された予ひずみ区分に対応する材料特性値と、今回のシェル要素の番号ｊとを互いに関連付けるデータがメモリ２２に格納される。したがって、前記スプリングバック解析プログラムは、そのデータを参照することにより、各シェル要素ごとに選択された材料特性値を用いることにより、材料のスプリングバックを解析することができる。
【０２７４】
続いて、Ｓ２０４において、今回の初期ひずみε^Ｐ _０（ｊ）に対して必要な修正が行われる。具体的には、上記Ｓ２０３において選択された材料特性値に厳密に対応する予ひずみε^Ｐ _ｇが厳密に今回の初期ひずみε^Ｐ _０（ｊ）に一致するとは限らないという事情を考慮し、一致しない場合には、一致するように今回の初期ひずみε^Ｐ _０（ｊ）が修正される。さらに、今回のシェル要素の番号ｊと、修正された初期ひずみε^Ｐ _０（ｊ）とを互いに関連付けるデータがメモリ２２に格納される。
【０２７５】
したがって、前記スプリングバック解析プログラムは、そのデータを参照することにより、各シェル要素ごとに修正された材料特性値を用いることにより、材料のスプリングバックを解析することができる。
【０２７６】
その後、Ｓ２０５において、番号ｊが最大値ｊＭＡＸ以上であるか否か、すなわち、材料の複数のシェル要素についてすべて材料特性値が選択されたか否かが判定される。今回は、番号ｊが最大値ｊＭＡＸ以上ではないと仮定すれば、判定がＮＯとなり、Ｓ２０６において、番号ｊが１増加させられ、その後、Ｓ２０２に移行する。以下、Ｓ２０２ないしＳ２０６が、新たなシェル要素について実行される。
【０２７７】
Ｓ２０２ないしＳ２０６の実行が何回か繰り返された結果、Ｓ２０５の判定がＹＥＳとなれば、この材料特性値決定プログラムの一回の実行が終了する。
【０２７８】
次に、前記バックストレス計算プログラムの内容を詳述する。
【０２７９】
図４０には、このバックストレス計算プログラムがフローチャートで示されている。まず、Ｓ３０１において、材料の前記複数のシェル要素に順に付与された番号ｋが１に設定される。次に、Ｓ３０２において、メモリ２２から、今回のシェル要素ｋに対応する初期応力−ひずみ関係が読み出される。この初期応力−ひずみ関係は、材料の今回のシェル要素ｋにプレス成形直後に残存する応力とひずみとの関係を意味しており、前記成形シミュレーション・プログラムによりシミュレートされたものである。
【０２８０】
なお、材料を仮想的に複数のシェル要素に分割する仕方は、バックストレスの計算と前記材料特性値の決定とで異ならないが、説明の便宜上、材料特性値の決定においては各シェル要素の番号がｊで表されているに対して、バックストレスαの計算においてはｋで表されている。
【０２８１】
その後、Ｓ３０３において、その読み出された初期応力−ひずみ関係に基づき、今回のシェル要素ｋのバックストレスαが計算される。さらに、その計算されたバックストレスαと今回のシェル要素の番号ｋとを互いに関連付けるデータがメモリ２２に格納される。したがって、前記スプリングバック解析プログラムは、そのデータを参照することにより、前記材料特性値のみならず、各シェル要素ごとに計算されたバックストレスαをも用いることにより、材料のスプリングバックを解析することができる。
【０２８２】
ここで、バックストレスαの計算手法を説明する。
【０２８３】
プレス成形シミュレーションにおいては、プレス成形されるべき材料が仮想的に複数のシェル要素に分割されるとともに、各シェル要素について平面応力状態が仮定される。したがって、各シェル要素においては、図４１において式（１０１）で示すように、板厚方向のすべての応力成分が無視される。この式（１０１）において「σ」は垂直応力を意味し、「τ」はせん断応力を意味している。
【０２８４】
このように仮定された各シェル要素においては、偏差応力ベクトル（σ’ｘｘ，σ’ｙｙ，σ’ｘｙ）の方向と、求めるべきバックストレスαの偏差成分（α’ｘｘ，α’ｙｙ，α’ｘｙ）の方向とが同じであると仮定すると、式（１０２）が成立する。この式（１２）において「Ａ」は定数を意味する。
【０２８５】
一方、偏差応力ベクトルσ’ｉｊについては、式（１０３）が成立する。この式（１０３）において「σｉｊ」は応力（垂直応力σのみならずせん断応力τも表す）、「δｉｊ」はクロネッカーのデルタ、「σｍ　」は静水圧下における応力σｉｊをそれぞれ意味する。応力σｍ　は、式（１０４）により求めることができる。
【０２８６】
さらに、バックストレスαに関しては、式（１０５）が成立する。したがって、上記定数Ａが、図４２の式（１０６）により求められる。この式（１０６）において「ε」は、垂直ひずみεとせん断ひずみγとを選択的に表す。具体的には、式（１０６）の右辺における「σ」が垂直応力σを表す場合には垂直ひずみεを表し、これに対し、せん断応力τを表す場合にはせん断ひずみγを表す。
【０２８７】
一方、上記式（１０２）により、図４２の式（１０７）で表される関係が成立する。さらに、バックストレスαの偏差成分α’については、偏差応力ベクトルσ’と同様に、式（１０８）が成立する。この式（１０８）において「αｉｊ」はバックストレス、「δｉｊ」はクロネッカーのデルタ、「αｍ　」は静水圧下におけるバックストレスαｉｊをそれぞれ意味する。応力αｍ　は、式（１０９）により求めることができる。
【０２８８】
したがって、それら式（１０６）ないし式（１０９）を用いることにより、バックストレスαが計算される。図４３には、求めるべきバックストレスαの相当値が応力σの相当値との関係においてグラフで示されている。
【０２８９】
以上のようにしてバックストレスαが計算された後、Ｓ３０４において、番号ｋが最大値ｋＭＡＸ以上であるか否か、すなわち、材料の複数の要素についてすべてバックストレスαが計算されたか否かが判定される。今回は、番号ｋが最大値ｋＭＡＸ以上ではないと仮定すれば、判定がＮＯとなり、Ｓ３０５において、番号ｋが１増加させられ、その後、Ｓ３０２に移行する。以下、Ｓ３０２ないしＳ３０５が、新たな要素について実行される。
【０２９０】
Ｓ３０２ないしＳ３０５の実行が何回か繰り返された結果、Ｓ３０４の判定がＹＥＳとなれば、このバックストレス計算プログラムの一回の実行が終了する。
【０２９１】
図４４には、本実施形態のシェル要素モデルを用いて成形シミュレーションを行ってスプリングバックを解析することによる効果がグラフで表されている。このグラフは、図４５で示すように定義されるスプリングバックの計算値と、そのための計算時間とが示されている。図４４には、スプリングバックの実験値が示されている。
【０２９２】
図４４には、さらに、比較のため、ソリッド要素モデルを用いた例と、従来例のシェル要素モデルを用いた例とが示されている。本実施形態のシェル要素モデルを用いれば、従来例のシェル要素モデルを用いた場合とほぼ同等な計算時間で、ソリッド要素モデルを用いた場合の計算精度に十分に匹敵する計算精度が実現される。
【０２９３】
以上の説明から明らかなように、本実施形態においては、図２のステップＳ４が前記（１）項における「成形シミュレーション工程」の一例を構成しているのである。
【０２９４】
さらに、本実施形態においては、図１３および図１４のＳ５０１ないしＳ５１４が前記（２）または（６）項における「暫定値解析工程」の一例を構成し、Ｓ５１５ないしＳ５１８が同項における「最終値解析工程」の一例を構成しているのである。
【０２９５】
さらに、本実施形態においては、図１４のＳ５１５ないしＳ５１８が前記（３）項における「最終値決定工程」の一例を構成しているのである。
【０２９６】
さらに、本実施形態においては、図１４のＳ５１８が前記（４）または（８）項における「工程」の一例を構成しているのである。
【０２９７】
さらに、本実施形態においては、図１４のＳ５１６が前記（５）または（９）項における「工程」の一例を構成しているのである。
【０２９８】
さらに、本実施形態においては、図１４のＳ５１６が前記（７）項における「応力近似関数特定工程」の一例を構成し、Ｓ５１７が同項における「仮面内応力決定工程」の一例を構成し、Ｓ５１５およびＳ５１８が互いに共同して同項における「最終値決定工程」の一例を構成しているのである。
【０２９９】
さらに、本実施形態においては、図１３のＳ５０２およびＳ５０３ならびに図２０のひずみ近似関数特定ルーチンが互いに共同して前記（１２）項における「ひずみ解析工程」の一例を構成しているのである。
【０３００】
さらに、本実施形態においては、図２０のＳ６０３が前記（１６）項における「ひずみ近似関数特定工程」の一例を構成し、図１３および図１４のＳ５０２ないしＳ５１８が同項における「面内応力解析工程」の一例を構成しているのである。
【０３０１】
さらに、本実施形態においては、図２のステップＳ７が前記（１７）項における「スプリングバック解析工程」の一例を構成しているのである。
【０３０２】
さらに、本実施形態においては、図２のステップＳ１が前記（１８）項における「実験値取得工程」の一例を構成し、ステップＳ２が同項における「材料特性値取得工程」の一例を構成しているのである。
【０３０３】
以上、本発明の一実施形態を図面に基づいて詳細に説明したが、これは例示であり、前記〔課題を解決するための手段および発明の効果〕の項に記載された態様を始めとして、当業者の知識に基づいて種々の変形，改良を施した形態で本発明を実施することが可能である。
【図面の簡単な説明】
【図１】本発明の一実施形態であるスプリングバック解析方法を実施するのに好適なスプリングバック解析装置を示す系統図である。
【図２】上記スプリングバック解析方法を示す工程図である。
【図３】板材をＵ字状に曲げることが可能なプレス成形機の一例を示す正面断面図である。
【図４】上記スプリングバック解析方法において応力−ひずみ関係の実験値が予ひずみごとに取得される様子を説明するためのグラフである。
【図５】図２におけるステップＳ４のスプリングバック解析に使用される材料硬化モデルおよびそのスプリングバック解析に必要な材料特性値を説明するためのグラフである。
【図６】上記スプリングバック解析に必要な材料特性値を説明するための別のグラフである。
【図７】上記スプリングバック解析に使用される材料硬化モデルを等方硬化モデルおよび実際の応力−ひずみ関係と対比しつつ説明するためのグラフである。
【図８】図２におけるステップＳ２の材料特性値取得を行うためにコンピュータにより実行される材料特性値取得プログラムの内容を概念的に表すフローチャートである。
【図９】図８の材料特性値取得プログラムにより材料特性値が取得される理論および原理を説明するためのいくつかの式を示す図である。
【図１０】図８におけるＳ１０６の詳細を材料軟化係数算出ルーチンとして概念的に表すフローチャートである。
【図１１】図８の材料特性値取得プログラムにより予ひずみと材料特性値との関係が離散的に取得される様子を表形式で示す図である。
【図１２】図２におけるステップＳ４の成形シミュレーションを行うためにコンピュータにより実行される成形シミュレーション・プログラムの内容を概念的に表すフローチャートである。
【図１３】図１２におけるＳ４０３の詳細を要素応力計算ルーチンとして概念的に表すフローチャートの一部である。
【図１４】上記要素応力計算ルーチンを表すフローチャートの別の一部である。
【図１５】図１３および図１４の要素応力計算ルーチンに利用される理論および原理を説明するためのいくつかの式を示す図である。
【図１６】図１３および図１４の要素応力計算ルーチンに利用される理論および原理を説明するためのグラフである。
【図１７】従来例のシェル要素に利用される変位−ひずみ関係式を示す図である。
【図１８】図１７の変位−ひずみ関係式により表現されるひずみ分布を板材の曲げ変形と共に示す図である。
【図１９】断面回転角θｙと板材の曲げ内面からの距離ｄｓとの関係が曲げの厳しさＬｂに応じて変化する様子を示すグラフである。
【図２０】図２におけるステップＳ４の成形シミュレーションを実行するための成形シミュレーション・プログラムに組み込まれたひずみ近似関数特定ルーチンの内容を概念的に表すフローチャートである。
【図２１】断面回転角θｙと板材の曲げ内面からの距離ｄｓとの関係についてソリッド要素モデルを用いて計算された値の精度を実測値と対比して検証するためのグラフである。
【図２２】図２０のひずみ近似関数特定ルーチンの実行結果のいくつかの例を示すグラフである。
【図２３】断面回転角θｙと板材の曲げ内面からの距離ｄｓとの関係についての３種類の計算値であってソリッド要素モデルを用いた計算値と上記実施形態のシェル要素モデルを用いた計算値と従来例のシェル要素モデルを用いた計算値とから成るものを示すグラフである。
【図２４】上記実施形態のシェル要素を用いて計算されたせん断ひずみγｘｚの分布と従来例のシェル要素を用いて計算されたせん断ひずみγｘｚの分布とを板材の曲げ変形と共に示す図である。
【図２５】上記実施形態のシェル要素モデルに利用される変位−ひずみ関係式を示す図である。
【図２６】上記実施形態のシェル要素を用いて計算された垂直ひずみεｘおよびせん断ひずみγｘｚの板厚方向における分布を板材の曲げ変形と共に示す図である。
【図２７】上記実施形態のシェル要素を用いて計算された垂直応力σｚの分布と従来例のシェル要素を用いて計算された垂直応力σｚの分布とを板材の曲げ変形と共に示す図である。
【図２８】図１４のＳ５１６において応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式の一部を示す図である。
【図２９】上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式の別の一部を示す図である。
【図３０】上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。
【図３１】上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。
【図３２】上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。
【図３３】上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。
【図３４】上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。
【図３５】垂直応力σｚを考慮しない垂直応力σｘ，σｚと垂直応力σｚを考慮した垂直応力σｘ’，σｙ’との関係を説明するためのいくつかの式を示す図である。
【図３６】図１４におけるＳ５１８において中立面を移動させるために利用される理論および原理を説明するためのグラフである。
【図３７】垂直応力σｚを考慮しない場合と考慮した場合と垂直応力σｚのみならず断面力のつり合いも考慮した場合とで垂直応力σｘが異なる様子を説明するためのグラフである。
【図３８】曲げ方向応力σｘと板材の曲げ内面からの距離ｄｓとの関係をソリッド要素モデルを用いて計算した場合と上記実施形態のシェル要素を用いて計算した場合と従来例のシェル要素を用いて計算した場合とについてそれぞれ示すグラフである。
【図３９】図２におけるステップＳ５の材料特性値決定を行うためにコンピュータにより実行される材料特性値決定プログラムの内容を概念的に表すフローチャートである。
【図４０】図２におけるステップＳ６のバックストレス計算を行うためにコンピュータにより実行されるバックストレス計算プログラムを概念的に表すフローチャートである。
【図４１】図４０のバックストレス計算プログラムによりバックストレスを計算するために利用される理論および原理を説明するためのいくつかの式の一部を示す図である。
【図４２】上記バックストレスを計算するために利用される理論および原理を説明するためのいくつかの式の別の一部を示す図である。
【図４３】図３８のバックストレス計算プログラムの実行内容を説明するためのグラフである。
【図４４】上記実施形態のシェル要素モデルを用いて板材のスプリングバックを計算した場合のその計算値と計算時間とを、ソリッド要素モデルを用いた場合と、従来例のシェル要素モデルを用いた場合と比較して示すグラフである。
【図４５】図４４におけるスプリングバックの定義を説明するための正面断面図である。
【符号の説明】
１２　コンピュータ
１６　外部記憶装置
３０　ＣＤ−ＲＯＭ
３２　ＭＯ
１００　プレス成形機
１１０　ポンチ
１１２　ダイス
１１４　パッド

Claims (15)

1. 曲げ変形させられるべき板材の両面の一方である曲げ外面をダイスとパッドとのうちの少なくともダイスによって支持しつつ他方の面である曲げ内面にポンチを加圧状態で接触させることによって前記板材をプレス成形する際のその板材の変形をコンピュータによるシミュレーションによって解析する方法であって、
   前記板材がそれの面方向に並んだ複数のシェル要素の集合体として前記コンピュータ上でモデル化されたシェル要素モデルを用い、前記各シェル要素の板厚方向応力を考慮した上で各シェル要素の面方向の応力である面内応力を解析する成形シミュレーション工程を含む成形シミュレーション解析方法。
2. 前記成形シミュレーション工程が、
   前記各シェル要素につき、前記板厚方向応力が０である平面応力状態を想定して前記面内応力の暫定値を解析する暫定値解析工程と、
   前記各シェル要素につき、前記解析された暫定値に基づき、かつ、前記板厚方向における力のつり合いを考慮することにより、前記面内応力の最終値を解析する最終値解析工程と
   を含む請求項１に記載の成形シミュレーション解析方法。
3. 前記最終値解析工程が、前記面内応力の、前記各シェル要素の断面上における総和である断面力が、その断面に隣接したシェル要素についての断面力と実質的に一致するように前記最終値を決定する最終値決定工程を含む請求項２に記載の成形シミュレーション解析方法。
4. 前記最終値決定工程が、前記各シェル要素の断面上において、前記板厚方向応力を考慮した面内応力の総和である断面力が前記平面応力状態における面内応力の総和である断面力に実質的に維持されるように前記各シェル要素の中立面の位置を決定し、その決定された中立面位置に基づき、前記最終値を決定する工程を含む請求項３に記載の成形シミュレーション解析方法。
5. 前記最終値解析工程が、前記解析された暫定値と、前記プレス成形時に前記板材が前記曲げ内面において前記ポンチから受ける接触圧とに基づき、前記最終値を解析する工程を含む請求項２ないし４のいずれかに記載の成形シミュレーション解析方法。
6. 前記成形シミュレーション工程が、
   前記各シェル要素につき、前記板厚方向応力が０である平面応力状態を想定して前記面内応力の暫定値を解析する暫定値解析工程と、
   前記各シェル要素につき、前記解析された暫定値に基づき、前記板厚方向応力の前記板厚方向における分布を近似的に定義する応力近似関数を特定し、その特定された応力近似関数を用いて前記面内応力の最終値を解析する最終値解析工程と
   を含む請求項１に記載の成形シミュレーション解析方法。
7. 前記最終値解析工程が、
   前記解析された暫定値に基づき、前記応力近似関数を特定する応力近似関数特定工程と、
   その特定された応力近似関数を用い、前記板厚方向応力と共に成立する面内応力を仮面内応力として決定する仮面内応力決定工程と、
   その決定された仮面内応力に基づき、前記面内応力の、前記各シェル要素の断面上における総和である断面力が、その断面に隣接したシェル要素についての断面力と実質的に一致するように前記最終値を決定する最終値決定工程と
   を含む請求項６に記載の成形シミュレーション解析方法。
8. 前記最終値決定工程が、前記各シェル要素の断面上において、前記板厚方向応力を考慮した面内応力の総和である断面力が前記平面応力状態における面内応力の総和である断面力に実質的に維持されるように前記各シェル要素の中立面の位置を決定し、その決定された中立面位置に基づき、前記最終値を決定する工程を含む請求項７に記載の成形シミュレーション解析方法。
9. 前記最終値解析工程が、前記解析された暫定値と、前記プレス成形時に前記板材が前記曲げ内面において前記ポンチから受ける接触圧とに基づき、前記応力近似関数を特定する工程を含む請求項６ないし８のいずれかに記載の成形シミュレーション解析方法。
10. 前記応力近似関数が、前記各シェル要素につき、前記板厚方向における力のつり合いを前提として前記分布を定式化する請求項６ないし９のいずれかに記載の成形シミュレーション解析方法。
11. 前記各シェル要素が複数の節点により定義され、かつ、それら複数の節点が、対応するシェル要素を幾何学的に特定するのに不可欠ではない節点を含まない請求項１ないし１０のいずれかに記載の成形シミュレーション解析方法。
12. 前記成形シミュレーション工程が、前記プレス成形による曲げの厳しさに関して前記板材と等価な等価板材がそれの面方向と板厚方向との双方に並んだ複数のソリッド要素の集合体として前記コンピュータ上でモデル化されたソリッド要素モデルを用い、前記曲げ変形前に前記等価板材の長手方向軸線に垂直である平面断面上の各点が前記曲げ変形によって回転する角度である断面回転角を前記等価板材の板厚方向位置に依存するように解析し、その断面回転角の解析結果を用い、前記等価板材のひずみを解析するひずみ解析工程を含む請求項１ないし１１のいずれかに記載の成形シミュレーション解析方法。
13. 前記ひずみ解析工程が、前記等価板材につき、前記断面回転角の解析結果に基づき、前記ひずみの前記板厚方向における分布を近似的に定義するひずみ近似関数を特定するひずみ近似関数特定工程を含み、
    前記成形シミュレーション工程が、前記各シェル要素につき、前記シェル要素モデルを前記特定されたひずみ近似関数と共に用い、かつ、前記板厚方向応力を考慮した上で前記面内応力を解析する面内応力解析工程を含む請求項１２に記載の成形シミュレーション解析方法。
14. 請求項１ないし１３のいずれかに記載の成形シミュレーション解析方法を実施するためにコンピュータにより実行されるプログラム。
15. 請求項１４に記載のプログラムをコンピュータ読み取り可能に記録した記録媒体。

Images