JP2006155254A - 応力解析方法、プログラムおよび記録媒体 - Google Patents

Abstract

【課題】プレス成形による板材の曲げ変形をコンピュータによって解析するために、シェル要素モデルを用いて各シェル要素における応力を解析する技術において、その解析に必要な計算時間の短縮を図りつつその解析の精度を向上させる。
【解決手段】各シェル要素につき、板厚方向応力が０である平面応力状態を想定して面内応力の暫定値を算出し（Ｓ５１２）、その暫定値に基づき、板厚方向応力の板厚方向における分布を近似的に定義する関数を特定する（Ｓ５１５）。その関数を用いて、板厚方向応力と共に力学的に成立する仮面内応力を算出し（Ｓ５１６）、面内応力の、各シェル要素の断面上における総和である断面力が増加する量を断面力増分として算出する（Ｓ５１７）。その断面力増分を、板厚方向に並んだ複数の位置にそれぞれ配分する配分量を決定し（Ｓ５１８）、その配分量で仮面内応力を補正して面内応力の最終値を算出する（Ｓ５１９）。
【選択図】図１４

Description

本発明は、プレス成形による板材の曲げ変形をコンピュータによって解析する技術に関するものであり、特に、その板材を近似的に表現する複数のシェル要素の集合体であるシェル要素モデルを用いることにより、各シェル要素における応力を解析する技術に関するものである。

ポンチ、ダイス等、金型を用いて板材を塑性加工によって成形する技術としてプレス成形が存在する。このプレス成形は、種々の製品を製造するために用いられ、そのような製品の一例として自動車がある。この自動車の製造工程においては、車体のパネルやフレーム、ドアのパネル、サスペンションのアーム等、種々の部品を製造するためにプレス成形が行われる。

このプレス成形は、一般に、曲げ変形させられるべき板材の両面の一方である曲げ外面をダイスとパッドとのうちの少なくともダイスによって支持しつつ他方の面である曲げ内面にポンチを加圧状態で接触させることによって行われる。このプレス成形中、板材には、曲げによる応力と接触による応力と摩擦力による応力が発生する。

このプレス成形によって製造された製品の形状精度が不良となる原因として、プレス成形後の製品のスプリングバックがある。このスプリングバックという現象は、プレス成形等、塑性加工において発生し、材料が負荷された後に除荷されると、材料が弾性回復してその形状が変化する現象である。この場合に、部分的に塑性変形を伴うことがある。

このスプリングバックの量は、金型の形状を含むプレス成形条件の調整によって改善することが可能である。そのプレス成形条件には、例えば、ポンチおよびダイスの各肩Ｒや、それらポンチとダイスとの間の半径方向および幅方向のクリアランスが含まれる。

したがって、プレス成形による製品の形状精度を向上させるためには、そのスプリングバック量を精度よく予測して金型の設計を行うことが有効である。板材の板厚ｔが厚いという理由や、板材の曲げ半径Ｒが大きいという理由によってその板材のスプリングバック量が大きいほど、そのスプリングバック量を精度よく予測し、その結果を金型の設計にフィードバックすることが強く要望される。

そのスプリングバック量を精度よく予測するために、プレス成形によって板材が変形する様子をコンピュータを用いた数値シミュレーションによって解析することが既に行われている。その数値シミュレーションは成形シミュレーションとも言われる。

その成形シミュレーションは、板材を、それの面方向と板厚方向との双方に沿って複数のソリッド要素に仮想的に分割してモデル化し、そのソリッド要素モデルのもとにＦＥＭ等の数値解析を行う手法と、板材を、それの面方向に複数のシェル要素に仮想的に分割してモデル化し、そのシェル要素モデルのもとにＦＥＭ等の数値解析を行う手法とに分類される。

一般に、先の手法は、板材変形の予測精度を容易に向上させることが可能である点で後の手法より優れ、一方、後の手法は、数値解析に必要な計算時間を容易に短縮可能な点で先の手法より優れている。したがって、予測精度と計算時間との両方において実用的に優れた手法の開発が要望されている。

従来のシェル要素を用いた成形シミュレーション解析は、プレス成形中に板材に板厚方向に発生する応力を考慮することができないが、薄板成形のように、板材の曲げの厳しさ（板材の板厚ｔに対するその板材の曲げ半径Ｒの比率）が比較的弱い場合には、その板材の変形を精度よく計算することができる。しかし、本発明者らは、この従来のシェル要素を用いた成形シミュレーションでは厚板の変形を精度よく予測するのは困難であることに気が付いた。

そこで、本発明者らは、板材がそれの面方向に並んだ複数のシェル要素の集合体としてコンピュータ上でモデル化されたシェル要素モデルを用いることにより、各シェル要素の板厚方向応力を考慮した上で各シェル要素の面方向における応力である面内応力を解析し、その面内応力の解析結果を用いて板材の変形を解析する技術を提案した。この技術は特許文献１に開示されている。

この開示技術によれば、板材の変形がシェル要素モデルを用いて解析されるため、ソリッド要素モデルを用いる場合に比較し、板材の変形をシミュレーションによって解析するための計算時間を短縮することが容易となる。しかも、この開示技術によれば、前述の従来のシェル要素モデルを用いて板材の変形を解析する場合とは異なり、その解析が板厚方向応力を考慮して行われるため、その解析の精度を、ソリッド要素モデルを用いる場合に近づけることが容易となる。このように、この開示技術によれば、計算時間の短縮と解析精度の向上とを両立させることが容易となる。

以下、この開示技術の内容をさらに具体的に説明するが、その説明の便宜上、板材においてその中立面が延びる方向すなわち曲げ方向をｘ方向、板厚方向をｚ方向といい、それら両方向に直角な方向をｙ方向（断面方向ともいう）という。それらのうちｘ方向とｙ方向とはいずれも、板材にとっての面内方向であり、残りのｚ方向は、板材の面と交差するという意味において面外方向ということが可能である。

シェル要素モデルを利用する場合には、板材の板厚方向における応力である板厚方向応力σｚを直接に計算することは困難であるが、例えば部分円筒を成す微小要素を想定するなどして、曲げ変形させられるべき板材における板厚方向（半径方向）における力のつり合いを式で表現すれば、板厚方向応力σｚと、板材の曲げ方向における応力である曲げ方向応力σｘ（板材の断面方向における応力である断面方向応力σｙも同様）との関係を誘導することが可能である。その誘導された関係を利用すれば、シェル要素モデルを用いて計算された曲げ方向応力σｘ（断面方向応力σｙも同様）から板厚方向応力σｚを計算することが可能である。

このような知見に基づき、上述の開示技術においては、各シェル要素につき、板厚方向応力が０である平面応力状態を想定して面内応力の暫定値が解析され、各シェル要素につき、その解析された暫定値に基づき、板厚方向における力のつり合いを考慮することにより、面内応力の最終値が解析される。

したがって、この開示技術によれば、板厚方向応力を直接に計算することは困難であるシェル要素モデルを用いるにもかからわず、そのシェル要素モデルを用いて計算された面内応力を有効に用いることにより、板厚方向応力を考慮した面内応力を計算することが可能となる。

しかしながら、各シェル要素につき、面内応力の暫定値、すなわち、平面応力状態において板厚方向応力を無視して計算された面内応力とは異なる値を結果的に有するように、面内応力の最終値を決定すると、各シェル要素の断面上における面内応力の総和である断面力が、その断面に隣接したシェル要素についての断面力と異なってしまう可能性がある。このように、互いに隣接したシェル要素間において断面力に関する力のつり合いが乱されてしまう事態は、互いに隣接したシェル要素によって共有される節点における面内方向における力のつり合いに反する。

このような知見に基づき、本発明者らは、面内応力の、各シェル要素の断面上における総和である断面力が、その断面に隣接したシェル要素についての断面力と実質的に一致するように、面内応力の最終値を決定する技術を提案した。この技術も特許文献１に開示されている。

この開示技術によれば、平面応力状態で計算された面内応力が板厚方向応力を考慮して変更されるにもかかわらず、互いに隣接したシェル要素間において断面力の連続性を保持することが可能となる。

さらに、本発明者らは、各シェル要素の断面上において、板厚方向応力を考慮した面内応力の総和である断面力が平面応力状態における面内応力の総和である断面力に実質的に維持されるように各シェル要素の中立面の位置を決定し、その決定された中立面位置に基づき、面内応力の最終値を決定する技術も提案した。この技術も特許文献１に開示されている。

この開示技術によれば、各シェル要素の断面に作用する面内応力の総和である断面力が実質的に維持されるように中立面の位置が決定され、その結果、互いに隣接したシェル要素間において断面力の連続性を保持することが可能となる。

具体的には、この開示技術によれば、まず、板厚方向応力を考慮しないで面内応力の暫定値が算出され、次に、その算出された面内応力の暫定値を利用し、かつ、板厚方向応力を考慮して、面内応力の最終値が算出される。その算出された面内応力の最終値の、板材の各断面における総和である断面力（「軸力」ともいう。）の計算値は、板厚方向応力を考慮することに伴い、曲げ方向において増減する。ここに、「断面力の計算値が曲げ方向において増減する」とは、同じ板材内において互いに隣接した断面間において断面力（軸力）の計算値が互いに異なることを意味する。しかし、この開示技術によれば、その板材の弾塑性特性を定義するためにその板材に設定される中立面の位置を適当量移動させることにより、その断面力の計算値の増減が解消される。
特開２００４−４２０９８号公報

しかしながら、本発明者らは、この開示技術には改良の余地があることに気が付いた。以下、このことを具体的に説明する。

本発明者らは、この開示技術の一具体例として、中立面の移動量を短時間で計算することを可能にするために、板材における面内応力の板厚方向における分布が予測式によって固定的に定義される状態で、中立面の移動量を計算する第１の手法を開発した。

しかしながら、この第１の手法を採用する場合には、例えば、ハット曲げ成形等、大きな張力がかかりながら板材が曲げ変形させられるプレス成形において、中立面の移動によって断面力を精度よく補正するのに限界がある。そのため、この第１の手法を採用する場合には、面内応力の計算精度が低下し、ひいてはスプリングバック量の予測精度も低下してしまう可能性がある。

さらに、本発明者らは、上述の開示技術の別の具体例として、中立面の移動量を高精度で計算することを可能にするために、中立面の移動量を可変パラメータとして、隣接したシェル要素間における断面力の連続性が保持されるように、最適化計算のための収束計算によって中立面の移動量を計算したうえで、その面内応力の分布を計算する第２の手法を開発した。

この第２の手法を採用する場合には、中立面の移動量を精度よく算出して断面力を精度よく補正することが可能となるが、そのために必要な計算時間が増加してしまう可能性がある。

以上説明した事情を背景として、本発明は、プレス成形による板材の曲げ変形をコンピュータによって解析するために、その板材を近似的に表現する複数のシェル要素の集合体であるシェル要素モデルを用いることにより、各シェル要素における応力を解析する技術において、その解析に必要な計算時間の短縮を図りつつその解析の精度を向上させることを課題としてなされたものである。

本発明によって下記の各態様が得られる。各態様は、項に区分し、各項には番号を付し、必要に応じて他の項の番号を引用する形式で記載する。これは、本発明が採用し得る技術的特徴の一部およびそれの組合せの理解を容易にするためであり、本発明が採用し得る技術的特徴およびそれの組合せが以下の態様に限定されると解釈すべきではない。すなわち、下記の態様には記載されていないが本明細書には記載されている技術的特徴を本発明の技術的特徴として適宜抽出して採用することは妨げられないと解釈すべきなのである。

さらに、各項を他の項の番号を引用する形式で記載することが必ずしも、各項に記載の技術的特徴を他の項に記載の技術的特徴から分離させて独立させることを妨げることを意味するわけではなく、各項に記載の技術的特徴をその性質に応じて適宜独立させることが可能であると解釈すべきである。

（１） プレス成形による板材の曲げ変形をコンピュータによって解析するために、前記板材がそれの面方向に並んだ複数のシェル要素の集合体として前記コンピュータ上でモデル化されたシェル要素モデルを用いることにより、前記各シェル要素の板厚方向における応力である板厚方向応力を考慮した上で各シェル要素の面方向における応力である面内応力を解析する応力解析方法であって、
前記各シェル要素につき、前記板厚方向応力が０である平面応力状態を想定して前記面内応力の暫定値を算出する暫定値算出工程と、
その算出された面内応力の暫定値に基づき、前記板厚方向応力の前記板厚方向における分布を近似的に定義する応力近似関数を特定する応力近似関数特定工程と、
その特定された応力近似関数を用いることにより、前記板厚方向における各位置ごとに、前記板厚方向応力と共に力学的に成立する面内応力を、前記板厚方向応力を考慮した仮面内応力として算出する仮面内応力算出工程と、
前記板厚方向応力の前記板厚方向における分布に基づき、前記面内応力の算出に際して前記板厚方向応力を考慮することに起因して前記面内応力の、前記各シェル要素の断面上における総和である断面力が増加する量を断面力増分として算出する断面力増分算出工程と、
その算出された断面力増分を、前記板厚方向に並んだ複数の位置にそれぞれ配分する配分量を決定する配分量決定工程と、
前記板厚方向における各位置ごとに、前記算出された仮面内応力を、前記決定された配分量で補正することにより、前記面内応力の最終値を算出する最終値算出工程と
を含む応力解析方法。

この方法においては、各シェル要素につき、板厚方向応力が０である平面応力状態を想定して面内応力の暫定値が算出される。その算出された面内応力の暫定値に基づき、板厚方向応力の板厚方向における分布を近似的に定義する応力近似関数が特定される。その特定された応力近似関数を用いることにより、板厚方向における各位置ごとに、板厚方向応力と共に力学的に成立する面内応力が、板厚方向応力を考慮した仮面内応力として算出される。

この方法においては、さらに、板厚方向応力の板厚方向における分布に基づき、面内応力の算出に際して板厚方向応力を考慮することに起因して面内応力の各シェル要素の断面上における総和である断面力が増加する量が断面力増分として算出される。その断面力増分は、板厚方向応力を考慮することに伴って発生する。その算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ配分されるように、その断面力増分の配分量が決定される。

この方法においては、さらに、板厚方向における各位置ごとに、算出された仮面内応力が、決定された配分量で補正されることにより、面内応力の最終値が算出される。

この方法においては、板厚方向応力を考慮することに伴う断面力増分が互いに隣接したシェル要素間において非一様に発生しないように、断面力増分の配分量の決定およびその配分量を用いた仮面内応力の補正を行うことが可能である。

したがって、この方法によれば、互いに隣接したシェル要素間において面内応力に関する力のつり合いが、板厚方向応力を考慮する前後を通じて、実質的に成立するように、面内応力を解析することが容易となる。

この方法によれば、仮面内応力の補正を、前述の開示技術に従って板材の中立面を移動させることによって行う場合に比較し、簡単な計算で行うことが可能となり、その結果、面内応力の解析に必要な全時間を短縮することが容易となる。したがって、この方法によれば、面内応力の解析精度を確保しつつ、その解析のための計算時間の短縮を図ることが容易となる。

この方法の用途の一例は、面内応力の解析結果に基づき、プレス成形中における板材の変形（プレス成形後における板材の形状を含む。）をコンピュータを用いてシミュレーションによって解析することであるが、この方法はそれ以外の用途に適用することが可能である。そのシミュレーションの解析結果の用途の一例は、その解析結果に基づいて板材のスプリングバックを解析することであるが、その解析結果はそれ以外の用途に利用することが可能である。

そのスプリングバックの解析結果の用途の一例は、当該板材に対して行われるプレス成形の条件を最適化することであるが、その解析結果はそれ以外の用途に利用することが可能である。そのプレス成形の条件には、例えば、そのプレス成形に用いられる金型の形状が含まれる。

本項および下記の各項に係る方法は、板材が厚板である場合に有効であるのはもちろんであるが、厚板の応力解析のみならず薄板の応力解析にも適用することが可能である。

本項および下記の各項において「応力」という用語は、垂直応力とせん断応力との少なくとも一方を包含し、「面内応力」という用語は、曲げ方向応力と断面方向応力との少なくとも一方を包含し、「ひずみ」という用語は、垂直ひずみとせん断ひずみとの少なくとも一方を包含する。

（２） 前記最終値算出工程は、前記板厚方向における各位置ごとに、前記算出された仮面内応力から、前記決定された配分量を引き算することにより、前記算出された仮面内応力を補正する（１）項に記載の応力解析方法。

板厚方向における各位置ごとの面内応力については、塑性力学的には、板厚方向応力を考慮して算出された値は、考慮しないで算出された値より、同じ各位置に作用する板厚方向応力と同じ量増加する。一方、板材における曲げ方向または断面方向における力のつり合いは、その板材において互いに隣接した断面間に、各断面全体に作用する面内応力の総和すなわち断面力について成立する。

したがって、板厚方向応力を考慮する前後を通じて、断面力に関して力のつり合いが維持されるようにするためには、板厚方向応力を考慮した仮面内応力の、各断面における総和である断面力から、板厚方向応力の、同じ各断面における総和である断面力増分を引き算すればよい。

前述の開示技術に従って板材の中立面の位置を補正することは、板材の曲げ方向または断面方向における力のつり合いを、板厚方向応力を考慮する前後を通じて維持するために有効である。しかし、中立面の位置は板材の断面形状によって幾何学的に決まる中央面の位置に固定する一方で、断面力増分を用いて仮面内応力を補正することも、板材の曲げ方向または断面方向における力のつり合いを、板厚方向応力を考慮する前後を通じて維持するために有効である。

ただし、その断面力増分は、断面全体についての面内応力の増分であって、板厚方向における各位置ごとの面内応力の増分ではないため、その断面力増分が判明しただけでは、板厚方向における各位置ごとの面内応力を精度よく補正することはできない。

とはいえ、板厚方向応力の板厚方向における分布を、板材の曲げ方向または断面方向における力のつり合いが少なくとも断面力に関して維持されるという条件と、板厚方向応力が実際に板厚方向において示す分布にできる限り近似する分布を示すという条件とが一緒に成立するように、予測することが可能である。

さらに、その予測された板厚方向応力の分布に従い、断面力増分を板厚方向における複数の位置にそれぞれ配分する配分量を決定し、その決定された配分量を仮面内応力から引き算すれば、その仮面内応力が精度よく補正されることになる。

このような知見に基づき、本項に係る方法においては、板厚方向における各位置ごとに、算出された仮面内応力から、決定された配分量が引き算されることにより、算出された仮面内応力が補正される。

したがって、この方法によれば、仮面内応力から断面力増分の配分量を引き算するという比較的簡単な計算により、その仮面内応力を補正することが可能となる。

（３） 前記配分量決定工程は、前記算出された断面力増分が前記複数の位置にそれぞれ互いに実質的に均等に配分されるように前記配分量を決定する（１）または（２）項に記載の応力解析方法。

前記（１）または（２）項に係る応力解析方法は、算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ、互いに均等ではないように配分されるように、断面力増分の配分量を決定する態様で実施することが可能である。この態様においては、例えば、板厚方向における位置に応じて線形的に変化するように配分量を決定したり、非線形的に変化するように配分量を決定することが可能である。

これに対し、本項に係る方法においては、算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ互いに実質的に均等に配分されるように、断面力増分の配分量が決定される。

したがって、この方法によれば、算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ互いに均等ではないように配分される場合に比較し、断面力増分の配分量を決定するために必要な計算を単純化してその計算に必要な時間を短縮することが容易となる。

ところで、例えば、板材の曲げ半径Ｒ［ｍｍ］をその板材の板厚ｔ［ｍｍ］で割り算して得られる曲げ指標が３より小さいというように、板材の曲げの厳しさが強い条件で実施されるプレス成形を想定すると、そのプレス成形により、当該板材は、それの板厚方向におけるほぼ全域において塑性域に到達する。一方、ひずみと板厚方向応力との関係を表すグラフの勾配を塑性域と弾性域とについて互いに比較すると、塑性域における方が弾性域におけるより緩やかである。このことは、板厚方向位置に対する板厚方向応力の依存性が、塑性域において弾性域におけるより小さいことを意味する。

したがって、上述のように板材の曲げの厳しさが強い条件で実施されるプレス成形による板材を解析対象として、それの面内応力を解析することが必要である場合には、板厚方向応力の板厚方向における分布を一様分布で近似することが妥当である。よって、前述の断面力増分の配分量も、板厚方向において一様であるように決定することが妥当である。

以上説明した知見に基づき、本項に係る方法においては、算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ互いに実質的に均等に配分されるように、断面力増分の配分量が決定される。

（４） 前記配分量決定工程は、前記板厚方向における位置に依存しないように前記配分量を決定する（１）ないし（３）項のいずれかに記載の応力解析方法。

前記（１）または（２）項に係る応力解析方法は、算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ、板厚方向における位置に依存するように、すなわち、位置に応じて異なるように配分されるように、断面力増分の配分量を決定する態様で実施することが可能である。これに対し、本項に係る方法においては、算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ、板厚方向における位置に依存しないように配分されるように、断面力増分の配分量が決定される。

したがって、この方法によれば、算出された断面力増分が、板厚方向に並んだ複数の位置にそれぞれ、板厚方向における位置に依存するように配分される場合に比較し、断面力増分の配分量を決定するために必要な計算を単純化してその計算に必要な時間を短縮することが容易となる。

（５） 前記断面力増分算出工程は、前記応力分布関数を前記板厚方向における位置に関して実質的に積分することにより、前記断面力増分を算出する（１）ないし（４）項のいずれかに記載の応力解析方法。

前記（１）ないし（４）項のいずれかに係る方法は、断面力増分を算出するために、板厚方向応力の板厚方向における分布を、前記（１）項における応力近似関数を用いることなく、定義し、そのようにして定義された分布に基づき、断面力増分を算出する態様で実施することが可能である。

これに対し、本項に係る方法においては、前記応力近似関数が、前記仮面内応力を算出するためのみならず、断面力増分を算出するためにも利用される。すなわち、この方法においては、その応力分布関数が前記板厚方向における位置に関して実質的に積分されることにより、断面力増分が算出されるのである。この断面力増分の算出は、例えば、各シェル要素ごとに行われる。

したがって、この方法によれば、仮面内応力と断面力増分とをそれらに共通の関数を用いて算出することが可能となるため、互いに異なる関数を用いて算出する場合に比較し、その算出のためのプログラム設計等、設計負担が軽減される。

（６） 前記プレス成形は、前記板材の曲げ半径Ｒ［ｍｍ］をその板材の板厚ｔ［ｍｍ］で割り算して得られる曲げ指標が３より小さい条件で実施される（１）ないし（５）項のいずれかに記載の応力解析方法。

本発明者らは、実験により、前記（１）ないし（５）項のいずれかに係る方法は、板材の曲げ半径Ｒ［ｍｍ］をその板材の板厚ｔ［ｍｍ］で割り算して得られる曲げ指標が３より小さい条件で実施されるプレス成形中に板材に発生する応力であっても高精度で解析することができることを確認した。

この実験結果に基づき、本項に係る方法は、上記曲げ指標が３より小さい条件で実施されるプレス成形中に板材に発生する応力を解析するために実施される。

（７） 前記プレス成形は、厚板成形とヘミング加工とハット曲げ加工との少なくとも一つを含む（１）ないし（６）項のいずれかに記載の応力解析方法。

本項における「厚板成形」は、例えば、板厚が２［ｍｍ］程度以上である板材を成形するために行われる。また、「ヘミング加工」は、例えば、自動車のドアのインナパネルとアウタパネルとを縁同士において互いに結合するために行われ、一般に、アウタパネルのうちインナパネルの縁から突き出た部分をそのインナパネルを挟むようにカールするために行われる。また、「ハット曲げ加工」は、深絞り成形の一例である。

（８） （１）ないし（７）項のいずれかに記載の方法を実施するためにコンピュータにより実行されるプログラム。

このプログラムがコンピュータにより実行されれば、前記（１）ないし（７）項のいずれかに係る方法と同じ作用効果が実現され得る。

本項における「プログラム」は、それの機能を果たすためにコンピュータにより実行される指令の組合せのみならず、各指令により処理されるファイルやデータをも含むように解釈することが可能である。

また、本項における「プログラム」は、各項に係る機能を果たすためにあえて作成された専用のプログラムとして、または専用のサブプログラム（モジュール等）の組合せとして作成することは不可欠ではなく、専用のサブプログラムと汎用のサブプログラムとの組合せとして作成することが可能であり、場合によっては、専用のサブプログラムを伴わない複数の汎用のサブプログラムの組合せとして作成する場合もあり得る。

（９） （８）項に記載のプログラムをコンピュータ読み取り可能に記録した記録媒体。

この記録媒体に記録されているプログラムがコンピュータにより実行されれば、前記（１）ないし（７）項のいずれかに係る方法と同じ作用効果が実現され得る。

本項における「記録媒体」は種々の形式を採用可能であり、例えば、フレキシブル・ディスク等の磁気記録媒体、ＣＤ、ＣＤ−ＲＯＭ等の光記録媒体、ＭＯ等の光磁気記録媒体、ＲＯＭ等のアンリムーバブル・ストレージ等の少なくとも一つを採用可能である。

以下、本発明のさらに具体的な実施の形態のうちの一つを図面に基づいて詳細に説明する。

図１には、本発明の一実施形態に従う応力解析方法を含むスプリングバック解析方法を実施するのに好適なスプリングバック解析装置（以下、単に「解析装置」という。）が示されている。

この解析装置は、入力装置１０とコンピュータ１２と出力装置１４と外部記憶装置１６とを備えている。入力装置１０はマウス，キーボード等を含むように構成される。コンピュータ１２はＣＰＵ等、プロセッサ２０と、ＲＯＭ，ＲＡＭ，ハードディスク等、メモリ２２と、それらプロセッサ２０とメモリ２２とを互いに接続するバス２４とを含むように構成される。出力装置１４はディスプレイ，プリンタ，プロッタ等として構成される。

外部記憶装置１６は、ＣＤ−ＲＯＭ３０，書き込み可能なＭＯ（光磁気ディスク）３２等、記録媒体が装填可能となっていて、装填状態においては、記録媒体に対するデータの読み取りおよび書き込みが必要に応じて行われる。

本実施形態においては、スプリングバックの解析に必要なプログラムがＣＤ−ＲＯＭ３０に記憶され、また、その解析に必要なデータがＭＯ３２に記憶されている。スプリングバックの解析に必要なプログラムには、材料特性値取得プログラムと、成形シミュレーション・プログラムと、材料特性値決定プログラムと、バックストレス計算プログラムと、スプリングバック解析プログラムとがある。

スプリングバックの解析時には、それらＣＤ−ＲＯＭ３０およびＭＯ３２から必要なプログラムおよびデータが読み出されてコンピュータ１２のＲＡＭまたはハードディスクに転送され、その後、プロセッサ２０によりそのプログラムが実行される。

図２には、このスプリングバック解析方法が工程図で表されている。

概略的に説明すれば、このスプリングバック解析方法は、板材としての厚板鋼板をプレス成形した後にその板材に生じるスプリングバックを解析するために実施される。そのプレス成形は、塑性ひずみが０．１５，０．２０というように、比較的高い塑性ひずみが生じる部分が存在するように行われる。

なお付言すれば、以下、「応力」という用語は、特に断りがない限り、垂直応力とせん断応力との双方を意味し、同様に、「ひずみ」という用語も、垂直ひずみとせん断ひずみとの双方を意味する。

図３には、プレス成形機１００が例示されている。このプレス成形機１００は、ワークＷとしての厚板をしわ押さえを利用せずにＵ字状に曲げ加工することが可能なものである。このプレス成形機１００においては、上ラム１０２と下ラム１０４とがすきまを隔てて互いに対向させられ、それら上ラム１０２と下ラム１０４との間において、ポンチ１１０とダイス１１２とパッド１１４とがワークＷを挟むように配置される。

ポンチ１１０は、ダイス１１２に接近する向きの力を上ラム１０２から受けて、予め設定された速度でワークＷに押し付けられる。ダイス１１２はステージ１２０に固定されている。パッド１１４は、油圧シリンダ１２２を介して下ラム１０４に支持されていて、ポンチ１１０に接近する向きの力を油圧シリンダ１２２から受けるようにされている。ダイス１１２はダイセット１２４によってステージ１２０に取り付けられる。

このスプリングバック解析方法においては、プレス成形されるべき材料と材質が実質的に同じ試験片を用いることにより、その材料の応力−ひずみ関係の実験値が複数の離散値として取得され、その取得された複数の実験値に基づき、その材料のスプリングバックを解析するのに必要な材料特性値が取得される。この材料特性値は、材料のバウシンガ効果が精度よく表現されるようにその材料の応力−ひずみ関係を定義するのに必要な情報であるということもできる。

また、このスプリングバック解析方法においては、プレス成形直後に材料に残存する応力−ひずみ関係を解析するために、有限要素法を用いることにより、材料のプレス成形による変形をシミュレートする成形シミュレーションが前記成形シミュレーション・プログラムの実行によって行われる。

このスプリングバック解析方法においては、その成形シミュレーション・プログラムの実行終了後、前記スプリングバック解析プログラムの実行により、材料のスプリングバックが解析される。この解析においては、有限要素法が用いられるとともに、

（１）前記材料特性値取得プログラムの実行により取得された材料特性値と、
（２）前記成形シミュレーション・プログラムの実行により計算された残留応力および残留ひずみ、板材の各要素ごとの中立面位置、ならびに板材の各要素ごとの板厚と、
（３）前記バックストレス計算プログラムの実行により計算されたバックストレスと

が用いられる。

さらに具体的に説明すれば、本実施形態においては、図２に示すように、まず、ステップＳ１において、作業者により、プレス成形されるべき材料と材質が実質的に同じ試験片が準備され、その後、図示しない試験機を用いることにより、その試験片に対して一軸応力状態において引張−圧縮試験が行われる。

この試験は具体的には、試験片を引張方向に負荷して塑性変形させた後に除荷し、さらに、逆方向すなわち圧縮方向に負荷して塑性変形させるというものであり、その間にその試験片の応力−ひずみ関係の実験値が複数の離散値として取得される。

それら複数の実験値は、試験片の除荷開始時における塑性ひずみである予ひずみを変化させるごとに取得される。図４には、除荷および圧縮時における応力−ひずみ関係を表すグラフが５つの予ひずみについてそれぞれ示されている。

次に、図２におけるステップＳ２において、コンピュータ１２により、前記材料特性値取得プログラムが実行され、それにより、上記取得された複数の実験値に基づいて材料特性値が各予ひずみごとに決定される。材料特性値取得プログラムの内容については、後に詳述する。

その後、ステップＳ３において、作業者により、前記成形シミュレーションを実行するのに必要なデータ、すなわち、数値解析用データが作成される。数値解析用データは、プレス成形すべき板材を仮想的に複数のシェル要素に分割するためのデータと、プレス成形時に板材に付与される荷重（板材にポンチ１１０およびパッド１１４から負荷される荷重を含む。）に関する条件と、複数のシェル要素に関する境界条件とを含んでいる。

本実施形態においては、各シェル要素を定義する複数の節点が、各シェル要素を幾何学的に特定するのに不可欠ではない節点は含まないようにされている。このようなシェル要素の一例は、四辺形を成し、それの４つの頂点を４つの節点とする四辺形要素である。各シェル要素を上記のように定義することにより、板材に設定される節点の数の増加が回避され、ひいては、節点数の増加に伴う計算時間の増加も回避される。

続いて、ステップＳ４において、コンピュータ１２により、上記成形シミュレーション・プログラムが実行され、それにより、プレス成形直後に板材に残存する応力−ひずみ関係が解析される。さらに、プレス成形直後における材料の形状も解析される。

この成形シミュレーション・プログラムは、概略的に説明すれば、材料の応力−ひずみ関係を近似する材料硬化モデルとして等方硬化モデルを使用する形式である。この成形シミュレーション・プログラムは、板材の各シェル要素の番号ｊに関連付けて、プレス成形直後に各シェル要素に存在する塑性ひずみである初期ひずみを残留ひずみとしてメモリ２２に格納する。この成形シミュレーション・プログラムの詳細については、後述する。

この成形シミュレーション・プログラムは、ＬＳＴＣ社により製造され、「ＬＳ−ＤＹＮＡ３Ｄ」という名称で日本総合研究所により販売されたプログラム、およびＥＳＩ社により製造され、「ＰＡＭ−ＳＴＡＭＰ」という名称で販売されたプログラムを用いるとともに、それらに新たなステップが追加されて構成されている。

その後、ステップＳ５において、コンピュータ１２により、前記材料特性値決定プログラムが実行され、それにより、上記ステップＳ２において取得された複数の材料特性値を用いることにより、板材の各シェル要素の成形に伴う塑性ひずみに対応する材料特性値が決定される。材料特性値決定プログラムの内容については、後に詳述する。

続いて、ステップＳ６において、コンピュータ１２により、前記バックストレス計算プログラムが実行され、それにより、板材の各シェル要素ごとにバックストレスαが計算される。バックストレス計算プログラムの内容については、後に詳述する。

その後、ステップＳ７において、コンピュータ１２により、前記スプリングバック解析プログラムが実行され、それにより、スプリングバックが解析される。

このスプリングバック解析プログラムは、前述の説明から明らかなように、

（１）前記材料特性値決定プログラムの実行により板材の各シェル要素ごとに決定された材料特性値と、
（２）前記成形シミュレーション・プログラムの実行により各シェル要素ごとに計算された残留応力および残留ひずみ、中立面位置ｒ_０ならびに板厚ｔ_０と、
（３）前記バックストレス計算プログラムの実行により各シェル要素ごとに計算されたバックストレスと

を用いることにより、板材のスプリングバックを解析するためのプログラムである。

このスプリングバック解析プログラムの一例は、日本総合研究所により製造されて「ＪＯＨ−ＮＩＫＥ３Ｄ」という名称で販売されたものである。

以上で、図２に示すスプリングバック解析方法の一回の実行が終了する。

前記スプリングバック解析プログラムは、スプリングバックを解析すべき材料の応力−ひずみ関係を種々の方式で定義可能に設計されている。その方式の一つに、応力−ひずみ関係を、平行四辺形を用いた線形の材料硬化モデルで近似する方式がある。図５には、試験片を引張方向に負荷（正負荷）して塑性変形させた後に除荷し、引き続いて逆向きすなわち圧縮方向に負荷（逆負荷）して再度塑性変形させた場合に、塑性ひずみε^Ｐと応力σとが変化する様子が平行四辺形で示されている。

ここで、図５における各種記号を説明すれば、「Ｙ_０」は、材料が初めての塑性変形を開始するときの降伏応力である初期降伏応力を示している。「Ｙ_１」は、材料が再降伏するときの応力である再降伏応力を示している。「Ｈ」は、塑性硬化係数を示している。塑性硬化係数Ｈは、正負荷時および逆負荷時の塑性域における応力−ひずみ関係を表す直線グラフ（後述の第３直線部５４）の勾配を角度φで表した場合にｔａｎφで表される。「β」は、再降伏応力Ｙ_１が初期降伏応力Ｙ_０ に対して低下する程度を規定する材料軟化係数を示している。「ε^Ｐ _０ 」は、除荷開始時における塑性ひずみε^Ｐである初期ひずみを示している。

前記スプリングバック解析プログラムは、応力−ひずみ関係を近似する材料硬化モデルとして、等方硬化モデルと移動硬化モデルと複合硬化モデルとの中から適当なものを選択できる。

この選択は、上記材料軟化係数βの値を決定することによって行われる。具体的には、図５に示すように、β＝０のときには、移動硬化モデルが選択され、このとき、再降伏応力Ｙ_１は、

Ｙ_１＝−Ｙ_０＋Ｈ・ε^Ｐ _０

なる式で定義される。また、β＝１のときには、等方硬化モデルが選択され、このとき、再降伏応力Ｙ₁ は、

Ｙ_１＝−Ｙ_０−Ｈ・ε^Ｐ _０

なる式で定義される。また、０＜β＜１のときには、複合硬化モデルが選択され、このとき、再降伏応力Ｙ₁ は、

Ｙ_１＝−Ｙ_０＋Ｈ・（１−２β）・ε^Ｐ _０

なる式で定義される。

図６のグラフは、プレス成形されるべき材料に関するものであるが、試験片に関しても同じグラフが描かれることになる。そして、試験片の応力−ひずみ関係を表すグラフを平行四辺形で近似する場合には、試験片の、正負荷時の弾性域における応力−ひずみ関係を表す第１直線部５０と、除荷時または逆負荷時の弾性域における応力−ひずみ関係を表す第２直線部５２とが互いに平行となるとともに、試験片の、正負荷時の塑性域における応力−ひずみ関係を表す第３直線部５４と、逆負荷時の塑性域における応力−ひずみ関係を表す第４直線部５６とが互いに平行となるように応力−ひずみ関係が近似される。

一方、応力−ひずみ関係を表すグラフのうちバウシンガ効果が現れる領域は、除荷開始点Ｐ_Ｓ から再降伏点Ｐ_Ｙ までの弾性域と、その再降伏点Ｐ_Ｙ から逆負荷終了点までの塑性域とである。

したがって、本実施形態においては、それら２つの領域における応力−ひずみ関係を近似的に表す折れ線、すなわち、第２直線部５２と第４直線部５６とが互いに連結されたものの幾何学的特徴が取得されるとともに、その取得された幾何学的特徴に基づき、直接的には試験片、間接的には材料につき、初期降伏応力Ｙ_０と、再降伏応力Ｙ_１と、塑性硬化係数Ｈと、材料軟化係数βとが取得される。

具体的には、再降伏応力Ｙ_１は、上記折れ線の折れ点における応力σとして取得される。塑性硬化係数Ｈは、第３直線部５４が第４直線部５６と平行であることを利用することにより、第４直線部５６の勾配を角度φで取得してそれのｔａｎφを算出することによって取得される。材料軟化係数βは、初期降伏応力Ｙ_０と、再降伏応力Ｙ_１と、塑性硬化係数Ｈと、予ひずみε^Ｐ _ｇとを、前述の式に類似した式、すなわち、

Ｙ_１＝−Ｙ_０＋Ｈ・（１−２β）・ε^Ｐ _ｇ

なる式に代入してβについて解くことにより取得される。なお、初期降伏応力Ｙ_０は、正負荷時における実験値に基づいて取得される。

このスプリングバック解析プログラムは、それら初期降伏応力Ｙ_０，再降伏応力Ｙ_１，塑性硬化係数Ｈおよび材料軟化係数βの他に、材料の弾性係数Ｅ（ヤング率Ｅとせん断弾性係数Ｇとを選択的に表す。）を利用することにより、スプリングバックを解析するように設計されている。すなわち、このスプリングバック解析プログラムにおいては、それら初期降伏応力Ｙ_０，再降伏応力Ｙ_１，塑性硬化係数Ｈ，材料軟化係数βおよび弾性係数Ｅがそれぞれ「材料特性値」を構成しているのである。

ところで、図５のグラフの横軸に取られているのは、弾性ひずみと塑性ひずみとの和である合計ひずみεではなく、塑性ひずみε^Ｐである。したがって、同図のグラフは、例えば、正負荷時の弾性域において、応力σを表す軸に平行に延びている。これに対して、図６のグラフの横軸に取られているのは、合計ひずみεである。したがって、同図のグラフは、例えば、正負荷時の弾性域において、応力σを表す軸に対して傾斜している。よって、それら２つのグラフにおいて、第１ないし第４直線部５０，５２，５４，５６の対応関係はそれら２つの図に示すものとなる。

したがって、図６における第１直線部５０の勾配を角度γで表した場合にｔａｎγとして求めた値は、材料のヤング率Ｅを表すことになる。一方、第１直線部５０と第２直線部５２とが互いに平行であると仮定されている。したがって、本実施形態においては、第２直線部５２の勾配からヤング率Ｅが取得される。

その取得されたヤング率Ｅから、予め定められた規則に従い、同じ材料についてのせん断弾性係数Ｇが誘導される。これにより、同じ材料につき、弾性係数Ｅと総称されるヤング率Ｅとせん断弾性係数Ｇとが取得されることになる。

応力−ひずみ関係は、材料の初期ひずみε^Ｐ _０（試験片の予ひずみε^Ｐ _ｇに相当する。）によって変化する。応力−ひずみ関係を近似的に表す平行四辺形の幾何学的特徴が初期ひずみε^Ｐ _０によって変化するのであり、よって、上記材料特性値も初期ひずみε^Ｐ _０によって変化する。したがって、本実施形態においては、前述のように、試験片の予ひずみε^Ｐ _ｇを変化させるごとに引張−圧縮試験が繰り返される。

図７には、本実施形態で使用される線形の材料硬化モデルが実線グラフで示され、等方硬化モデルが二点鎖線グラフで示され、実際の応力−ひずみ関係が破線グラフで示されている。その実際の応力−ひずみ関係は、正負荷時の弾性域および塑性域においては等方硬化モデルで表される関係（二点鎖線）と一致し、逆負荷時の弾性域および初期の塑性域においては本実施形態の材料硬化モデルで表される関係（実線）と一致する。本実施形態で使用される材料硬化モデルは、線形であるが、予ひずみε^Ｐ _ｇに応じて形状が変化するように定義されている。なお、本実施形態で使用される線形の移動硬化モデルを表す平行四辺形は、同図においては、それの２つの鋭角を含む一対の部分が省略されている。

本実施形態においては、材料硬化モデルとして、バウシンガ効果を表現し得る移動硬化モデルまたは複合硬化モデルを使用することができるとともに、材料の各シェル要素について予想される初期ひずみε^Ｐ _０に応じて形状が異なる材料硬化モデルが使用される。ただし、各シェル要素について予想される初期ひずみε^Ｐ _０が異なっても、材料硬化モデルが平行四辺形で表現されることは維持される。

したがって、図７に示すように、本実施形態で使用される線形の材料硬化モデルは、実際の応力−ひずみ関係を表すグラフと全体において一致するわけではないが、バウシンガ効果が現れる領域、すなわち、除荷後に弾性変形する領域と塑性変形する領域とについては、十分に高い精度で一致する。

このような理由から、本実施形態においては、材料の応力−ひずみ関係が、平行四辺形を用いた線形の材料硬化モデルで近似させられているのであり、よって、本実施形態によれば、低ひずみ領域のみならず高ひずみ領域においてもバウシンガ効果を精度よくシミュレート可能となり、ひいては、スプリングバックを精度よく解析可能となる。

図７に示すように、前記折れ線は、除荷および逆負荷時における実際の応力−ひずみ関係を表すグラフにできる限り一致するように取得される。具体的には、そのグラフ上に、同一直線上に位置せず、かつ、相互に一定距離以上離れた３点を暫定的に設定することを、そのグラフのうち折れ線として抽出される可能性のある部分の一方の端点から他方の端点に向かって順次行い、それら暫定的な３点により構成される暫定的な折れ線の角度が極大となったときに、それら３点を最終的な３点に決定するとともに、それら最終的な３点により構成される折れ線を最終的な折れ線に決定する。

それら３点は、除荷開始点に最も近い第１実測点６０と、次に近い第２実測点６２と、最も遠い第３実測点６４とから構成される。それら３つの実測点６０，６２，６４における各応力−ひずみ関係は、実際の応力−ひずみ関係を表すグラフ上に位置することから、実測値と一致する。

なお、図７に示す例においては、折れ線を構成する２本の直線部のうち、除荷開始点（図示しない）に近い方は、前記平行四辺形の第２直線部５２と完全にではなく部分的に一致し、また、他方の直線は、平行四辺形の第４直線部５６と完全にではなく部分的に一致している。

ここで、前記材料特性値取得プログラムの内容を詳述する。

図８には、この材料特性値取得プログラムがフローチャートで表されている。まず、ステップＳ１０１（以下、単に「Ｓ１０１」で表す。他のステップについても同じとする。）において、実験値を表すデータがＭＯ３２からコンピュータ１２のメモリ２２に取り込まれる。

次に、Ｓ１０２において、予ひずみε^Ｐ _ｇに順に付与すべき番号ｉが１に設定される。その後、Ｓ１０３において、今回の予ひずみε^Ｐ _ｇ（ｉ）に関連する応力−ひずみ関係を表すデータがメモリ２２から取り込まれる。

続いて、Ｓ１０４において、その取り込まれたデータに基づき、今回の予ひずみε^Ｐ _ｇ（ｉ）に関連する応力−ひずみ関係を表すグラフが想定されるとともに、そのグラフが折れ線で近似させられる。それにより、平行四辺形上の前記３つの実測点６０，６２，６４が前述のようにして抽出される。その後、Ｓ１０５において、第２実測点６２と第３実測点６４とを結ぶ第２直線５２の勾配を角度φで算出することにより、今回の予ひずみε^Ｐ _ｇ（ｉ）に対応する塑性硬化係数Ｈ（ｉ）（＝ｔａｎφ）が算出される。

続いて、Ｓ１０６において、図９において式（１）で示す前述の式を用いることにより、今回の予ひずみε^Ｐ _ｇ（ｉ）に対応する材料軟化係数β（ｉ）が算出される。本実施形態においては、材料軟化係数βが０＜β＜１の範囲内となること、すなわち、材料硬化モデルが複合硬化モデルとなることを前提に算出される。以下、その算出方法を具体的に説明する。

上記式（１）を用いて材料軟化係数βを算出する際、図９において式（２）で示す第１条件と、すべての予ひずみε^Ｐ _ｇ について、式（３）で示す条件が成立するという第２条件とが制約条件として採用される。また、上記式（１）における「Ｙ_１」と「Ｈ」は、実験値により取得されたものを採用することにする。

なお、再降伏応力Ｙ_１については、種々の実験により、予ひずみε^Ｐ _ｇに関して単調増加を示すことが確認されており、予ひずみε^Ｐ _ｇ＝０のときには、

Ｙ_１＝−Ｙ_０

という関係が成立し、また、すべての予ひずみε^Ｐ _ｇについて、

Ｙ_１＞−Ｙ_０

という関係が成立する。

上記式（２）の第１条件を考慮し、かつ、上記式（１）を用いると、式（４）が得られる。この式（４）は式（５）に変形できる。この式（５）を用い、かつ、上記式（３）で示す第２条件を考慮すると、式（６）が得られ、この式（６）から式（７）が得られる。

また、上記式（３）で示す第２条件を考慮し、かつ、上記式（１）を「β」について解くと、式（８）が得られる。

そして、本実施形態においては、上記式（７）で示す条件が成立するか否かを判定し、成立しない場合には、成立するように今回の予ひずみε^Ｐ _ｇ（ｉ）が修正される。この修正は例えば、今回の予ひずみε^Ｐ _ｇ（ｉ）を、

１．０１×（−Ｙ_１／Ｈ）

として計算したり、実験値のうち上記（７）式で示す条件を満たすもののうち最も小さいものを選択することによって行われる。

上記式（７）で示す条件が成立したならば、必要に応じて修正された今回の予ひずみε^Ｐ _ｇ（ｉ）と、実験値により取得された初期降伏応力Ｙ_０，再降伏応力Ｙ_１および塑性硬化係数Ｈ（ｉ）とを上記式（１）に代入することにより、今回の材料軟化係数β（ｉ）を算出する。

さらに、その算出された材料軟化係数β（ｉ）が、前記式（２）で示す第１条件が成立するか否か、すなわち、０＜β＜１の範囲内にあるか否かを判定し、その第１条件が成立する場合には、その算出された材料軟化係数β（ｉ）が最終値とされる。その第１条件が成立しない場合には、その第１条件が成立するように材料軟化係数β（ｉ）が修正され、さらに、それに合わせて初期降伏応力Ｙ_０も修正される。

材料軟化係数β（ｉ）および初期降伏応力Ｙ_０の修正は例えば、もとの材料軟化係数β（ｉ）が０以下である場合には、材料軟化係数β（ｉ）を０より大きい設定値、例えば、０．１に変更するとともに、材料軟化係数β（ｉ）が０．１である場合の初期降伏応力Ｙ_０を上記式（１）を用いて計算する。その計算された初期降伏応力Ｙ_０が上記式（３）で示す条件を満たす場合には、材料軟化係数βの最新値β（ｉ）および初期降伏応力Ｙ_０がそれぞれそのまま最終値として採用される。

これに対して、変更された材料軟化係数β（ｉ）の下での初期降伏応力Ｙ_０が上記式（３）で示す条件を満たさない場合には、材料軟化係数β（ｉ）が例えば０．１増加させられ、同様にして、新たな初期降伏応力Ｙ_０が計算される。このような処理は、材料軟化係数βが１を超えない範囲内で繰り返され、その結果、材料軟化係数β（ｉ）が、実験値に基づく再降伏応力Ｙ_１および塑性硬化係数Ｈ（ｉ）の下に、上記式（１）ないし式（３）のすべてを満たすように修正されることになる。

以上、もとの材料軟化係数β（ｉ）が０以下である場合の修正方法を説明したが、１以上である場合にも、同様な方法で修正が行われる。

図１０には、Ｓ１０６の詳細が材料軟化係数算出ルーチンとしてフローチャートで表されている。

まず、Ｓ２５１において、前記抽出された第２実測点６２を用いることにより、再降伏応力Ｙ_１が決定される。次に、Ｓ２５２において、その決定された再降伏応力Ｙ_１と、前記算出された塑性硬化係数Ｈ（ｉ）と、今回の予ひずみε^Ｐ _ｇ（ｉ）とが図９の式（７）で示す条件を満たすか否かが判定される。

今回は、前記式（７）で示す条件が満たされると仮定すれば、Ｓ２５２の判定がＹＥＳとなり、Ｓ２５３において、上記決定された再降伏応力Ｙ_１と、前記算出された塑性硬化係数Ｈ（ｉ）と、実験値により取得された初期降伏応力Ｙ_０とを前記式（１）に代入することにより、今回の材料軟化係数β（ｉ）が算出される。その後、Ｓ２５４に移行する。

これに対して、今回は、前記式（７）で示す条件が満たされないと仮定すれば、Ｓ２５２の判定がＮＯとなり、Ｓ２５５において、今回の予ひずみε^Ｐ _ｇ（ｉ）が、前述のようにして、前記式（７）で示す条件を満たすように修正される。

その後、Ｓ２５３において、修正された予ひずみε^Ｐ _ｇ（ｉ）の下で前記式（１）を用いることにより今回の材料軟化係数β（ｉ）が算出される。続いて、Ｓ２５４に移行する。

いずれの場合にも、その後、Ｓ２５４において、その算出された材料軟化係数β（ｉ）が、図９の式（２）で示す条件を満たすか否かが判定される。今回は、前記式（２）で示す条件が満たされると仮定すれば、その判定がＹＥＳとなり、上記算出された材料軟化係数β（ｉ）と、実験値に基づく初期降伏応力Ｙ_０とがそのまま最終値として採用される。以上で本ルーチンの一回の実行が終了する。

これに対して、今回は、前記式（２）で示す条件が満たされないと仮定すれば、Ｓ２５４の判定がＮＯとなり、Ｓ２５６において、今回の材料軟化係数β（ｉ）が前述のようにして修正される。その後、Ｓ２５７において、修正された材料軟化係数β（ｉ）の下に、前記式（１）を用いることにより、初期降伏応力Ｙ_０が算出される。

続いて、Ｓ２５８において、その算出された初期降伏応力Ｙ_０が０より大きいか否かが判定される。今回は、０より大きいと仮定すれば、判定がＹＥＳとなり、その修正された材料軟化係数β（ｉ）と、その算出された初期降伏応力Ｙ_０とがそれぞれ最終値として採用される。以上で、本ルーチンの一回の実行が終了する。

また、今回は、その算出された初期降伏応力Ｙ_０が０以下であると仮定すれば、Ｓ２５８の判定がＮＯとなり、Ｓ２５６に戻る。このＳ２５６においては、材料軟化係数β（ｉ）が再度、前述のようにして修正され、その後、Ｓ２５７において、再度修正された材料軟化係数β（ｉ）の下に新たな初期降伏応力Ｙ_０が算出される。

続いて、Ｓ２５８において、前記の場合と同様にして、その算出された初期降伏応力Ｙ_０が０より大きいか否かが判定される。Ｓ２５６ないしＳ２５８の実行が何回か繰り返された結果、Ｓ２５８の判定がＹＥＳとなれば、最終的に修正された材料軟化係数β（ｉ）と初期降伏応力Ｙ_０とがそれぞれ最終値として採用される。以上で本ルーチンの一回の実行が終了する。

以上のようにしてＳ１０６の実行が終了すれば、その後、図８のＳ１０７において、第１実測点６０と第２実測点６２とをつなぐ第２直線５２の勾配が角度γで算出され、さらに、ｔａｎγとして、今回の予ひずみε^Ｐ _ｇ（ｉ）に対応するヤングＥ（ｉ）が算出される。その算出されたヤング率Ｅ（ｉ）から、予め定められた規則に従ってせん断弾性係数Ｇ（ｉ）が誘導される。このようにして、ヤング率Ｅ（ｉ）とせん断弾性係数Ｇ（ｉ）とが取得される。

続いて、Ｓ１０８において、番号ｉが最大値ｉＭＡＸ以上であるか否かが判定される。複数の実験値のすべてが材料特性値を決定するために利用されたか否かが判定されるのである。今回は、番号ｉが最大値ｉＭＡＸ以上ではないと仮定すれば、その判定がＮＯとなり、Ｓ１０９において、番号ｉが１増加させられ、その後、Ｓ１０３に移行する。以下、Ｓ１０３ないしＳ１０８が、新たな予ひずみε^Ｐ _ｇ（ｉ）に関連して実行されることになる。

Ｓ１０３ないしＳ１０９の実行が何回か繰り返された結果、Ｓ１０８の判定がＹＥＳとなれば、Ｓ１１０において、予ひずみε^Ｐ _ｇの変化領域が、設定された最小値と最大値との間において、一定間隔で複数の区分に分割される。

このとき、各区分に対応する予ひずみε^Ｐ _ｇに対応する材料特性値が存在しない場合には、材料特性値が存在しない予ひずみε^Ｐ _ｇに近接する複数の予ひずみε^Ｐ _ｇであって材料特性値が存在するものを内挿することにより、必要な材料特性値を生成する。その結果、予ひずみε^Ｐ _ｇと材料特性値との関係が、図１１に表形式で表されるように、予ひずみε^Ｐ _ｇの複数の代表値について離散的に取得されることになる。

以上で、この材料特性値取得プログラムの一回の実行が終了する。

次に、前記成形シミュレーション・プログラムの内容を詳述する。

図１２には、その成形シミュレーション・プログラムの内容がフローチャートで概念的に表されており、そのうちのステップＳ４０３の詳細が図１３および図１４に要素応力計算ルーチンとしてフローチャートで表されている。以下、この成形シミュレーション・プログラムの内容を説明するが、それに先立ち、その前提となる理論および原理を説明する。

この成形シミュレーション・プログラムにおいては、板材のプレス成形による変形が時間増分形の有限要素法を用いて解析される。

ここに、有限要素法は、よく知られているように、材料を複数の有限要素に分割し、各要素の各節点における応力と変位とを互いに関連付ける剛性マトリクスを用い、静的な問題に対しては各要素の各節点における外力と内力とのつり合いを条件に、連立方程式である剛性方程式の繰返し計算をそれの解が真の解に収束するまで行う数値解析法である。本実施形態においては、複数の要素剛性マトリクスが集合して成る全体剛性マトリクスの各要素が未知数であって、求めるべき解とされ、その解が真の解に収束するまで行う繰返し計算が行われる。

本実施形態においては、その有限要素法が増分計算法に組み合わされて実行される。増分計算法は、よく知られているように、非線形挙動や動的過程中の材料挙動を解析するために、時間を微小区間に分割し、それぞれの時間増分で計算された各物理量（変位、回転角、応力、ひずみ）の増分を順次累積して行くことによって、所定の時間までの解析を進める数値解析法である。

本実施形態においては、プレス成形中のポンチ１１０の一回の動作に伴う材料挙動が、各々微小時間を有する複数の計算ステップに分割されており、それら複数の計算ステップの終了により、一連の成形シミュレーションが終了する。さらに、本実施形態においては、増分計算法として動的な陽解法が採用されている。

この動的な陽解法により、ばねの運動方程式と等価な動的な運動方程式が解かれる。具体的には、各シェル要素の各節点の加速度が計算され、その加速度の時間積分によって各節点の速度が計算され、その速度の時間積分によって各節点の変位が計算される。

この成形シミュレーション・プログラムにおいては、さらに、弾塑性体のための応力更新アルゴリズムに基づき、応力およびひずみが解析される。この応力更新アルゴリズムは、フォン・ミーゼス(Von Mises)の降伏条件のもと、各シェル要素に時刻ｔから時刻ｔ＋Δｔまでの各回の計算ステップ中に生じたひずみ増分に対する応力増分が解析される。この応力更新アルゴリズムは、図１３および図１４を参照して後述される要素応力計算ルーチンとして具体化されている。

この応力更新アルゴリズムにおいては、今回の計算ステップにおける変位増分ｕと回転角増分θとが計算される。変位増分ｕは、本来、「ｄｕ」または「Δｕ」として表記すべきであるが、ここでは単に「ｕ」として表記する。また、変位増分ｕは、各成分の変位増分ｕ，ｖ，ｗの総称である。同様にして、回転角増分θは、本来、「ｄθ」または「Δθ」として表記すべきであるが、ここでは単に「θ」として表記する。また、回転角増分θは、各成分の回転角増分θｘ，θｙ，θｚの総称である。

具体的には、変位増分ｕの計算は、力のつり合いを考慮した増分型の全体剛性方程式を用いることによって行われる。これに対し、回転角増分θの計算は、モーメントのつり合いを考慮した増分型の全体剛性方程式を用いることによって行われる。いずれについても、全体剛性方程式は、未知の成分を有する全体剛性マトリクスを含み、かつ、その全体剛性マトリクスの未知の成分は、今回の計算ステップにおいては前回の計算ステップまでの計算値に等しくされる。本実施形態においては、全体剛性マトリクスの未知の成分が有限要素法によって求められる。

前記ひずみ増分は、等方ひずみ増分ε_ｍと偏差ひずみ増分ｅ_ｉｊとに分離される。その等方ひずみ増分ε_ｍは、各成分ごとのひずみ増分ε_ｉｊを用いて計算され、そのひずみ増分ε_ｉｊは、後述のひずみ近似関数を用いて計算される。

等方ひずみ増分ε_ｍから静水圧増分３Ｋε_ｍ（Ｋ：体積弾性係数）が求められる。この静水圧増分３Ｋε_ｍに前回の計算ステップまでの静水圧成分σ_ｉｊ ^０を加えることにより、静水圧成分σ_ｍが更新される。すなわち、更新された静水圧成分σ_ｍ ^ｔ＋Δｔが計算されるのであり、このことが図１５に式（１１）で示されている。

各成分ごとのひずみ増分ε_ｉｊから等方ひずみ増分ε_ｍを差し引くことにより、各成分ごとの偏差ひずみ増分ｅ_ｉｊが求められる。

この偏差ひずみ増分ｅ_ｉｊから偏差応力成分２Ｇｅ_ｉｊ（Ｇ：せん断弾性係数）が求められる。したがって、偏差応力成分２Ｇｅ_ｉｊは、平面応力状態において弾性域内において求められることになる。この偏差応力成分２Ｇｅ_ｉｊに前回の計算ステップまでの偏差応力成分Ｓ_ｉｊ ^０を加えることにより、試行応力（偏差応力Ｓの試行値）Ｓ_ｉｊ ^{Ｔｒｉａｌ}が求められる。すなわち、更新された試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}が計算されるのであり、このことが式（１２）で示されている。

この試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}を式（１３）に代入することにより、相当応力の試行解σ^{Ｔｒｉａｌ}が求められる。

この相当応力の試行解σ^{Ｔｒｉａｌ}が弾性域内にあり、各シェル要素の材料が未だ降伏していない場合には、上記試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}がそのまま、更新された偏差応力成分（今回の偏差応力成分）Ｓ_ｉｊ ^{Ｆｉｎａｌ}とされる。

一方、応力空間を表す図１６の（ａ）に示すように、相当応力の試行解σ^{Ｔｒｉａｌ}が塑性域に達すると、各シェル要素の材料が降伏したと判断される。この場合には、同図の（ｂ）に示すように、図１５に式（１４）で示すラジアルリターンアルゴリズム（半径引戻し解法）により、試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}が降伏曲面上へ引き戻される。

そのラジアルリターンアルゴリズムにおいては、相当塑性ひずみ増分ε^Ｐから、降伏曲面の半径を表す相当応力σと、最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とが計算される。さらに、それら両計算値が比較され、相当応力σに一致するときの最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値が最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の最終値に決定される。

なお付言すれば、前述のフォン・ミーゼスの降伏条件を採用する場合には、例えば、図１５の式（１６）を用いることより増分Δε^Ｐを計算し、その後、式（１４）を用いることにより最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}を計算することが可能である。

各シェル要素の材料が降伏した場合には、その引き戻された最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の最終値が、更新された偏差応力成分として求められる。

そして、図１５に式（１５）で示すように、上記更新された静水圧成分σ_ｍ ^ｔ＋Δｔと、上記更新された偏差応力成分Ｓ_ｉｊ ^{Ｆｉｎａｌ}とにより、各成分ごとに更新された応力σ_ｉｊ ^ｔ＋Δｔが計算される。

以上説明した応力更新の１回の計算ステップは、板材を構成するすべてのシェル要素について順に実行される。

従来のシェル要素では、曲げ変形の前後で断面の平面性が維持されると仮定されており、そのため、シェル要素の中立面は常に板厚中央面と一致するとともに、断面回転角θｙが、曲げ方向位置を表すｘ座標値には依存するが、板厚方向位置を表すｚ座標値には依存しない値とされる。

ここに、断面回転角θｙとは、シェル要素の中立面の法線に対するそのシェル要素の辺の角度を意味する。より正確には、曲げ変形前にシェル要素の長手方向軸線に垂直であった平面断面を変形前断面、曲げ変更後にその変形前断面が移動した断面を変形後断面とした場合に、その変形後断面上の任意の点上の接線と変形前断面との成す角度を意味する。

したがって、断面回転角θｙは、図１７の式（２１）で表され、曲げ方向の垂直ひずみεｘは式（２２）で表され、板厚方向のせん断ひずみγｘｚは式（２３）で表される。ただし、それらの式において、「ｕ」は、ｘ方向の変位を表し、「ｗ」は、ｚ方向の変位を表している。

よって、従来のシェル要素を用いる場合には、図１８に示すように、垂直ひずみεｘは板厚方向に線形に分布する一方、せん断ひずみγｘｚが板厚方向に沿って変化しないように分布する。

しかし、本発明者らの研究により、そのような仮定は、曲げ変形させられるべき板材が薄板である場合には妥当であるが、厚板である場合のように、板材の曲げの程度が厳しい場合には妥当性を失うことが判明した。

図１９には、断面回転角θｙと、板材内の、それの曲げ内面からの距離ｄｓとの関係が曲げの厳しさＬｂによって変化する様子がグラフで表されている。このグラフは、本発明者らの実験結果に基づいて作成されたものである。ここに、曲げの厳しさＬｂは、曲げ半径Ｒ［ｍｍ］を板厚ｔ［ｍｍ］で割り算した値である曲げ指標Ｒ／ｔを用い、その値が小さいほど曲げの厳しさが強くなるように表現される。

図１９のグラフによれば、曲げの厳しさＬｂが３以上である領域（曲げの厳しさが弱い領域であって、同図においては「Ｌｂ＝３」という表記に関連付けられた水平線の下側の領域）においては、板厚方向位置にかかわらず断面回転角θｙがほとんど変化せず、よって、従来のシェル要素についての仮定が妥当であることが分かる。しかし、曲げの厳しさＬｂが３より小さい領域（曲げの厳しさが強い領域であって、上記水平線の上側の領域）においては、板厚方向位置によって断面回転角θｙが顕著に変化し、よって、従来のシェル要素についての仮定が妥当ではないことが分かる。

そこで、本発明者らは、厚板中の任意の部位の断面回転角θｙをシミュレーションによって計算した。このシミュレーションにおいては、厚板が複数のソリッド要素に仮想的に分割される。ソリッド要素モデルを用いてシミュレーションが行われるのである。

このシミュレーションのために実行されるプログラムが前記成形シミュレーション・プログラムに組み込まれている。図２０には、そのプログラムがひずみ近似関数特定ルーチンとして概念的に表されている。

このひずみ近似関数特定ルーチンにおいては、まず、Ｓ６０１において、曲げ指標Ｒ／ｔが異なる複数の代表板材（潜在的な等価板材）のそれぞれにつき、有限要素法により、各ソリッド要素ごとに変位ｕ（各成分の変位ｕ，ｖ，ｗ）と回転角θ（各成分の回転角θｘ，θｙ，θｚ）とが計算される。次に、Ｓ６０２において、それら計算された変位ｕと回転角θとに基づき、かつ、図２５の式（３２）および式（３３）を用いることにより、ひずみがｚ座標値に関連付けて計算される。

続いて、Ｓ６０３において、その計算されたひずみとｚ座標値との関係に基づき、ひずみの板厚方向における分布を近似的に定義するひずみ近似関数が特定される。その後、Ｓ６０４において、その特定されたひずみ近似関数が、各板材の曲げ指標Ｒ／ｔに関連付けられる。以上で、このひずみ近似関数特定ルーチンの一回の実行が終了する。

このひずみ近似関数特定ルーチンは、スプリングバックを解析すべき板材ごとに実行することは不可欠ではない。このひずみ近似関数特定ルーチンの過去の実行により、実質的に同じ曲げの厳しさＬｂに関してひずみ近似関数が既に特定されている場合には、その特定されたひずみ近似関数を利用することにより、ひずみ近似関数特定ルーチンの新たな実行を省略することができるからである。

本発明者らは、ソリッド要素モデルを用いたシミュレーションの精度を検証するため、厚板につき、実験によって断面回転角θｙと曲げ内面からの距離ｄｓとの関係を測定した。さらに、本発明者らは、その測定された関係を表す実験値と上述のシミュレーションによる計算値とを互いに比較した。

図２１には、厚板につき、断面回転角θｙと距離ｄｓとの関係が実験値と計算値との双方によって示されている。同図から明らかなように、それら実験値と計算値との間には良好な一致が認められる。したがって、厚板につき、実験値に代えて計算値を用いることが可能であり、よって、実験なしでも、断面回転角θｙと距離ｄｓとの関係、すなわち、断面回転角θｙのｚ座標値への依存性を取得することが可能である。

図２２には、上記シミュレーションによる計算値から計算された板厚方向のせん断ひずみγｘｚ（前記ひずみの一例）と曲げ内面からの距離ｄｓとの関係が曲げ指標Ｒ／ｔの各値に関連付けてグラフで示されている。同図には、さらに、曲げ指標Ｒ／ｔの値がそれぞれ１．８３，１．６１，２．０５である場合のせん断ひずみγｘｚと距離ｄｓとの関係を近似するひずみ近似関数がグラフで示されている。このひずみ近似関数は、距離ｄｓを表す変数ｚを変数とする関数によってせん断ひずみγｘｚを定義する式である。このひずみ近似関数は、例えば、２次関数、３次関数、４次関数または５次関数である多項式として構成することが可能である。

このひずみ近似関数を利用すれば、シェル要素を用いた成形シミュレーションであっても、ソリッド要素を用いた場合とほぼ同等の精度で、ひずみの板厚方向における分布を計算することが可能となる。そして、図２３には、せん断ひずみγｘｚに関し、ソリッド要素を用いた成形シミュレーションによる計算値と、本実施形態のシェル要素を用いた成形シミュレーションによる計算値とがグラフで示されている。両者は良好な一致を示している。

図２３には、さらに、比較のために、従来例のシェル要素を用いた成形シミュレーションによる計算値も示されている。従来例のシェル要素においては、せん断ひずみγｘｚがｚ座標値に依存しない一定値であると仮定されている。

図２４には、本実施形態のシェル要素を用いて計算されたせん断ひずみγｘｚの分布と、従来例のシェル要素を用いて計算されたせん断ひずみγｘｚの分布とがグラフで表されている。本実施形態のシェル要素を用いて計算されたせん断ひずみγｘｚの分布は、前記ひずみ近似関数が２次関数である場合の一例を示している。このように、本実施形態のシェル要素を用いる場合には、従来例のシェル要素を用いる場合に比較し、厚板の材料特性の実際、すなわち、せん断ひずみγｘｚの分布の実際が精度よく再現される。

従来例のシェル要素とは異なり、本実施形態のシェル要素については、断面回転角θｙがｘ座標値のみならずｚ座標値にも依存するが、このことが図２５において式（３１）で表されている。

このような依存性により、本実施形態のシェル要素については、曲げ方向の垂直ひずみεｘが図２５の式（３２）で表され、板厚方向のせん断ひずみγｘｚは式（３３）で表される。

図２６には、本実施形態のシェル要素を用いた計算値が、せん断ひずみγｘｚと垂直ひずみεｘとのそれぞれの、板厚方向における分布に関して示されている。せん断ひずみγｘｚの分布は、前記ひずみ近似関数を２次より高次の関数である場合の一例を示している。本実施形態のシェル要素によれば断面回転角θｙをｚ座標値に依存させることが可能になり、その結果、せん断ひずみγｘｚのみならず垂直ひずみεｘについても実際の分布が精度よく再現される。

以上、本実施形態のシェル要素と従来例のシェル要素とをひずみに関して互いに比較したが、次に、板厚方向の垂直応力σｚに関して互いに比較する。

図２７に示すように、従来例のシェル要素においては、前述のように、垂直応力σｚが無視されて０であると仮定される。これに対し、本実施形態のシェル要素においては、垂直応力σｚが考慮される。

板厚方向の垂直応力σｚは、厚板の曲げ変形による寄与分σｚｂ（以下、「垂直応力σｚｂ」で表す。）と、厚板がポンチ１１０から受ける接触圧ｐによる寄与分σｚｃ（以下、「垂直応力σｚｃ」で表す。）との和であると考えることができる。垂直応力σｚｃは、接触圧ｐによって厚板に生ずるせん断変形に起因する。

垂直応力σｚｂは、厚板（厚肉円筒）内の微小要素のｚ方向（厚肉円筒の半径方向）における力のつり合いに着目することにより、垂直応力σｘから解析することが可能である。その力のつり合いは、極座標で記述するとともに時間的に離散化すれば、図２８における式（４１）で表される。この式（４１）における「ｒ」は、ｚ座標軸上の変数であり、「ｚ」に置き換えることが可能である。

本実施形態においては、まず、垂直応力σｚを無視して垂直応力σｘが暫定的に解析され、その暫定値に基づき、かつ、上記式（４１）を用いて、垂直応力σｚｂの板厚方向における分布を近似的に表す応力近似関数が特定される。垂直応力σｘの暫定値は、厚板の平面応力状態における垂直応力σｘに一致する。

垂直応力σｚｂのための応力近似関数の特定は具体的に次のようにして行うことが可能である。

まず、その応力近似関数の基本形を選択する。この基本形は、ｎ次関数により構成したり、指数関数により構成することが可能である。

図２８の式（４２）は、その応力近似関数の一例である。この式（４２）で表される応力近似関数において「ａ」は、厚板の曲げ内面のｚ方向位置を表し、また、「ｂ」は、厚板の曲げ外面のｚ方向位置を表している。中立面のｚ方向位置をｒ_０、板厚をｔ_０でそれぞれ表すと、「ａ」は式（４３）で、「ｂ」は式（４４）でそれぞれ表される。

したがって、上記式（４２）で表される垂直応力σｚｂは、厚板の両表面上において０となり、このことは、厚板の曲げ内面がポンチ１１０から受ける接触圧１１０が考慮されないことを示している。

垂直応力σｚｂを式（４２）で表される応力近似関数ｆ（ｚ）で表すこととすれば、前記式（４１）が式（４５）に変形される。

上記式（４２）は式（４６）で表すように展開できる。また、応力近似関数ｆ（ｚ）の一次導関数ｆ’（ｚ）は、式（４７）で表される。それら式（４６）および（４７）を上記式（４５）に代入すれば、式（４８）が得られる。

式（４８）中の各係数に式（４３）および式（４４）を代入するとともに、変数ａ０，ａ１，ｌ，ｍを用いた変数置換を行うと、式（４８）中の各係数が図２８の式（４９）で示すように整理される。この整理の結果を用いれば、図２８の式（４８）が図２９の式（５０）に変形される。

式（５１）で示す変数置換を行うと、式（５０）が式（５２）に変形される。この式（５２）においては、「Ｘ」と「Ｙ」と「Ｗ」と「Ｚ」とが既知である一方、「ｌ」と「ｍ」と「Ａ’」とが未知である。

ところで、本実施形態においては、各シェル要素が、それの板厚方向に積層された複数のレイヤに関連付けられるとともに、各シェル要素ごとに、垂直応力σｘの板厚方向における分布を表す関数（各レイヤの位置すなわちｚ座標値を変数とする）が予め定義される。それにより、本実施形態においては、各シェル要素につき、各ｚ座標値ごとに垂直応力σｘが求められ、それらｚ座標値と垂直応力σｘとの関係を前提にして最小２乗法が式（５２）に適用されることにより、未知数ｌ，ｍ，Ａ’が計算される。

図３０の式（５３）は、図２９の式（５２）について２乗和を表している。この式（５３）を未知数ｌ，ｍ，Ａ’に関してそれぞれ偏微分すれば、式（５４）が得られ、これは式（５５）に変形できる。この式（５５）で表される３つの連立方程式を解けば、未知数ｌ，ｍ，Ａ’が求められ、これにより、応力近似関数ｆ（ｚ）が特定される。

なお付言すれば、関数によって板厚方向分布が定義される変数に、垂直応力σｚを選んだり、せん断応力τｘｚを選んだり、両者を選ぶことが可能である。

以上、応力近似関数ｆ（ｚ）の例として４次関数である場合を説明したが、次に、３次関数である場合を説明する。

垂直応力σｚｂを図３１の式（５６）で表される応力近似関数ｆ（ｚ）で表すこととすれば、応力近似関数ｆ（ｚ）の一次導関数ｆ’（ｚ）は、式（５７）で表される。それら式（５６）および（５７）を前記式（４５）に代入すれば、式（５８）が得られる。この式（５８）を変数ｃに関して整理すれば、式（５９）が得られる。

式（６０）で示す変数置換を行うと、上記式（５９）が式（６１）に変形される。この式（６１）においては、「Ｘ」と「Ｙ」と「Ｚ」とが既知である一方、「ｃ」と「Ａ’」とが未知である。各シェル要素につき、各レイヤのｚ座標値と各レイヤの垂直応力σｘとの関係を前提にして最小２乗法が式（６１）に適用されることにより、未知数ｃ，Ａ’が計算される。式（６１）の２乗和が図３２の式（６２）で表されている。

式（６２）を未知数ｃに関して偏微分すれば、式（６３）が得られ、これは式（６４）に変形できる。これに対し、式（６２）を未知数Ａ’に関して偏微分すれば、式（６５）が得られ、これは式（６６）に変形できる。それら式（６４）および（６６）の連立方程式を解くことにより、未知数ｃ，Ａ’が求められ、これにより、応力近似関数ｆ（ｚ）が特定される。

図３３には、板厚方向の垂直応力σｚが、厚板の曲げ変形による垂直応力σｚｂと、接触圧ｐによる垂直応力σｚｃとの和に相当することが重み付け関数を用いて式（６９）で表されている。同図には、さらに、その垂直応力σｚｃの板厚方向における分布を近似的に表す応力近似関数が式（７０）で表されている。この応力近似関数を採用する場合には、式（７１）で表すように、厚板の曲げ内面（ｚ＝ａ）で垂直応力σｚｃが接触圧ｐに等しくなり、一方、曲げ外面（ｚ＝ｂ）で垂直応力σｚｃが０となる。

図３３には、式（６９）で表される重み付け関数の一例が式（７２）で表されている。この例においては、垂直応力σｚｂのための重み係数（式（６９）において「ｋ_１」に相当する。）と、垂直応力σｚｃのための重み係数（式（６９）において「ｋ_２」に相当する。）とがいずれも固定値として設定される。さらに、この例においては、それら重み係数およびが互いに等しく設定される。

図３３には、上記重み付け関数の別の例が式（７３）で表されている。この式（７３）における「α」が式（７４）で表されている。この例においては、垂直応力σｚｂのための重み係数（式（６９）において「ｋ_１」に相当する。）と、垂直応力σｚｃのための重み係数（式（６９）において「ｋ_２」に相当する。）とがいずれも、２次関数を用いて定義される可変値として設定される。

図３３には、上記重み付け関数のさらに別の例が式（７５）で表されている。この式（７５）における「ｘ」が式（７６）で表されている。この例においては、垂直応力σｚｂのための重み係数（式（６９）において「ｋ_１」に相当する。）と、垂直応力σｚｃのための重み係数（式（６９）において「ｋ_２」に相当する。）とがいずれも、飽和型関数（例えば、ブリルアン関数）を用いて定義される可変値として設定される。

本実施形態においては、板材に作用する垂直応力σｘ（＝第１の面内応力である曲げ方向応力）と垂直応力σｙ（＝第２の面内応力である断面方向応力）とがシェル要素を用いて解析されるにもかかわらずその解析が垂直応力σｚ（＝板厚方向応力）を考慮して行われる。ここに、垂直応力σｚは、図３３において式（６９）により、Ｆ（σｚｂ，σｚｃ））として表現される。

具体的には、本実施形態においては、まず、板厚方向応力σｚが０である平面応力状態が仮定されて曲げ方向応力σｘと断面方向応力σｙとが求められる。次に、それら求められた曲げ方向応力σｘと断面方向応力σｙ、または前回の計算ステップにおいて求められた曲げ方向応力σｘと断面方向応力σｙから、板厚方向応力σｚの板厚方向における分布を定義するために前記応力近似関数が同定される。その後、その同定された応力分布関数を用いることにより、板厚方向応力σｚが求められ、その求められた板厚方向応力σｚを前提にして、すなわち、板厚方向応力σｚを考慮して、曲げ方向応力σｘ’と断面方向応力σｙ’とが、前述のラジアルリターンアルゴリズムを用いて求められる。

プラントル−ルイス(Prandtl-Reuss)の流れ法則によれば、各成分の偏差応力Ｓについて図３４の式（８１）が成立する。また、垂直応力σｚを無視する場合には、各成分の偏差応力ｓについて式（８２）が成立する。これに対し、垂直応力σｚを考慮する場合には、各成分の偏差応力ｓについて式（８３）が成立する。式（８２）における「σｘ」および「σｙ」は、垂直応力σｚを考慮しない場合のｘ方向垂直応力（曲げ方向応力）およびｙ方向垂直応力（断面方向応力）をそれぞれ意味し、一方、式（８３）における「σｘ’」および「σｙ’」は、垂直応力σｚを考慮する場合のｘ方向垂直応力およびｙ方向垂直応力をそれぞれ意味している。

式（８２）中の第１式の右辺の値と式（８３）中の第１式の右辺の値とが互いに一致することから、式（８４）が誘導される。また、式（８２）中の第２式の右辺の値と式（８３）中の第２式の右辺の値とが互いに一致することから、式（８５）が誘導される。

それら式（８４）および式（８５）から、垂直応力σｘ’については式（８６）、垂直応力σｙ’については式（８７）がそれぞれ誘導される。それら式（８６）および式（８７）によれば、垂直応力σｚを考慮する垂直応力σｘ’，σｙ’は、考慮しない垂直応力σｘ，σｙに対して垂直応力σｚの分だけ増加することが分かる。この増加は、厚板の横断面全体に作用する垂直応力σｘ’，σｙ’の各総和である断面力（垂直応力σｚを考慮した断面力）が、垂直応力σｘ，σｙの各総和である断面力（平面応力状態での断面力）に対して垂直応力σｚの総和の分だけ増加することを意味する。

図３５には、プレス成形中に、板材を構成する各シェル要素に、現実には、各シェル要素の曲げと、各シェル要素の中立面に平行な軸力Ｆによる伸びとが重畳的に作用することが概念的に表されている。図３５には、さらに、そのような材料力学的挙動が、理論的には、各シェル要素の曲げと伸びとに分離されることも概念的に表されている。ここに、軸力Ｆは、各垂直応力σｘ，σｙごとに存在し、かつ、各垂直応力σｘ，σｙの、シェル要素の断面上における総和であることから、断面力とも称される。

本実施形態においては、垂直応力σｚを考慮した断面力Ｆが、平面応力状態での断面力に維持されるように、垂直応力σｘ’，σｙ’が、板厚方向における各位置ごとに補正される。すなわち、本実施形態においては、板材の中立面位置ｒ_０の変更を行うことなく、垂直応力σｘ’，σｙ’が補正され、それにより、垂直応力σｘ，σｙの最終値が垂直応力σｘ’’，σｙ’’として誘導される。

ところで、前述のように、板厚方向における各位置ごとの面内応力σｘ，σｙについては、塑性力学的には、板厚方向応力σｚを考慮して算出された値は、考慮しないで算出された値より、同じ各位置に作用する板厚方向応力σｚと同じ量増加する。一方、板材における曲げ方向および断面方向における力のつり合いは、その板材において互いに隣接した断面間に、各断面全体に作用する面内応力の総和すなわち断面力Ｆについて成立する。

したがって、板厚方向応力σｚを考慮する前後を通じて、断面力Ｆに関して力のつり合いが維持されるようにするためには、板厚方向応力σｚを考慮した仮面内応力σｘ’，σｙ’の、各断面における総和である断面力Ｆから、板厚方向応力σｚの、同じ各断面における総和である断面力増分Ｆｚを引き算すればよい。

その断面力増分Ｆｚは、図３６において式（９１）で表される。この式（９１）において「σｚ」は、前述の応力近似関数を意味する。

本実施形態においては、板材がそれの面方向に並んだ複数のシェル要素によって近似的に表現されており、各シェル要素は、板厚の中心のみを結ぶ面である。ただし、接触処理は、板厚を考慮して内外表面で行われる。図３７に示すように、要素の軸力とモーメントは、面内の積分点Ｉｐ１ないしＩｐ４（積分点の数は、使用要素により異なるが、典型的には４つである。）で計算され、それが形状関数を用いて節点が決定される。また、面内の各積分点は、板厚方向に何層かに積分点（図３７において「ｘ」で示す。）が存在する。各シェル要素の内部の連続領域を代表する離散的な複数の代表点が積分点として設定され、それら設定された積分点は、各シェル要素の内部において着目すべき評価点として利用することが可能である。

この場合には、それら積分点のうち、各シェル要素の板厚方向に並んだ複数の積分点のそれぞれについて垂直応力σｚを算出し、それら積分点について垂直応力σｚの算出値を合計することにより、上述の断面力増分Ｆｚが算出される。

断面力増分Ｆｚを用いて仮面内応力σｘ’，σｙ’を補正することは、板材の曲げ方向および断面方向における力のつり合いを、板厚方向応力σｚを考慮する前後を通じて維持するために有効である。ただし、その断面力増分Ｆｚは、断面全体についての面内応力σｘ，σｙの増分であって、板厚方向における各位置ごとの面内応力σｘ，σｙの増分ではないため、その断面力増分Ｆｚが判明しただけでは、板厚方向における各位置ごとの仮面内応力σｘ’，σｙ’を精度よく補正することはできない。

とはいえ、板厚方向応力σｚの板厚方向における分布を、板材の曲げ方向および断面方向における力のつり合いが少なくとも断面力Ｆに関して維持されるという条件と、板厚方向応力σｚが実際に板厚方向において示す分布にできる限り近似する分布を示すという条件とが一緒に成立するように、予測することが可能である。

さらに、その予測された板厚方向応力σｚの分布に従い、断面力増分Ｆｚを板厚方向における複数の位置にそれぞれ配分する配分量σｚｍを決定し、その決定された配分量σｚｍを仮面内応力σｘ’，σｙ’から引き算すれば、その仮面内応力σｘ’，σｙ’が精度よく補正されることになる。

ところで、例えば、板材の曲げ半径Ｒ［ｍｍ］をその板材の板厚ｔ［ｍｍ］で割り算して得られる曲げ指標が３より小さいというように、板材の曲げの厳しさＬｂが強い条件で実施されるプレス成形を想定すると、そのプレス成形により、当該板材は、それの板厚方向におけるほぼ全域において塑性域に到達する。一方、ひずみと板厚方向応力σｚとの関係を表すグラフの勾配を塑性域と弾性域とについて互いに比較すると、塑性域における方が弾性域におけるより緩やかである。このことは、板厚方向応力σｚの板厚方向位置に対する依存性が、塑性域において弾性域におけるより小さいことを意味する。

したがって、上述のように板材の曲げの厳しさＬｂが強い条件で実施されるプレス成形による板材を解析対象として、それの面内応力σｘ，σｙを解析することが必要である場合には、板厚方向応力σｚの板厚方向における分布を一様分布で近似することが妥当である。よって、前述の断面力増分Ｆｚの配分量σｚｍも、板厚方向において一様であるように決定することが妥当である。

このような知見に基づき、本実施形態においては、図３６において式（９２）で表すように、算出された断面力増分Ｆｚが、板厚ｔ_０で割り算されることにより、板厚方向において平均化され、その平均値が、断面力増分Ｆｚの配分量σｚｍとして決定される。具体的には、図３７に示す例においては、断面力増分Ｆｚが、板厚方向に並んだ複数の積分点の数によって割り算され、それにより、この断面力増分Ｆｚが、それら積分点にそれぞれ互いに均等に割り当てられることになる。

さらに、図３６において式（９３）で表すように、算出された配分量σｚが、板厚方向応力σｚを考慮した曲げ方向応力（仮面内応力）σｘ’から引き算されることにより、その曲げ方向応力σｘの最終値σｘ’’が誘導される。同様にして、式（９４）で表すように、算出された配分量σｚが、板厚方向応力σｚを考慮した断面方向応力（仮面内応力）σｙ’から引き算されることにより、その断面方向応力σｙの最終値σｙ’’が誘導される。

図３８には、曲げ方向応力σｘを例にとり、面内応力についての複数種類の計算値がそれぞれ複数のグラフで表されている。それら計算値は、

（１）シェル要素モデルを用い、かつ、板厚方向応力σｚが０であると仮定して計算された曲げ方向応力σｘと、
（２）シェル要素モデルを用い、かつ、板厚方向応力σｚを考慮して計算された曲げ方向応力σｘ’と、
（３）曲げ方向応力σｘ’から配分量σｚｍを引き算して計算された曲げ方向応力σｘの最終値σｘ’’（＝σｘ’−σｚｍ）と

を含んでいる。図３８には、さらに、板厚方向応力σｚの値もグラフで表されている。

図３８において、最終値σｘ’’が、曲げ方向応力σｘ’から配分量σｚｍを引き算して得られた値であるにもかかわらず、最終値σｘ’’を表すグラフが、曲げ方向応力σｘ’を表すグラフを応力の正方向に配分量σｚｍの絶対値と等しい量だけ平行移動させたグラフに相当するのは、この例においては、配分量σｚｍが圧縮応力であって、その符号が負であるからである。

したがって、本実施形態によれば、垂直応力σｚを考慮する場合と考慮しない場合とで、垂直応力σｘの総和である断面力が変化しない。これにより、同じ節点を共有する２つのシェル要素間において断面力が、垂直応力σｚを考慮したことが原因で不連続となってしまう不自然な事態が回避される。

以上、図１２に示す成形シミュレーション・プログラムの内容を概略的に説明したが、以下、具体的に説明する。

この成形シミュレーション・プログラムにおいては、まず、Ｓ４０１において、成形の解析に必要なデータがＭＯ３２から読み込まれるとともに、解析の初期化が行われる。

次に、Ｓ４０２において、板材に作用する荷重に関する境界条件が適用される。本実施形態においては、パッド１１４から板材に作用する荷重に関する境界条件が適用される。

続いて、Ｓ４０３において、板材を構成する複数のシェル要素のうち今回の実行対象である注目シェル要素に作用する応力が計算される。このＳ４０３の詳細が図１３および図１４に要素応力計算ルーチンとしてフローチャートで表されている。この要素応力計算ルーチンは、前記応力更新アルゴリズムを具体化するために設計されている。

この要素応力計算ルーチンにおいては、まず、Ｓ５０１において、前述のように、増分型の剛性方程式において全体剛性マトリクスをそれの各成分が前回値に等しい状態で用いることにより、各シェル要素の各節点における変位増分ｕおよび回転角増分θが計算される。これらの値をもとに、板厚中央面でのひずみが計算される。次に、Ｓ５０２において、前記ひずみ近似関数に基づき、各成分ごとにひずみ増分ε_ｉｊが計算される。ここに、ひずみ増分ε_ｉｊは、垂直ひずみとせん断ひずみとの双方を含んでいる。

その後、Ｓ５０３において、その計算されたひずみ増分ε_ｉｊに基づき、等方ひずみ増分ε_ｍが計算される。続いて、Ｓ５０４において、その計算された等方ひずみ増分ε_ｍに基づき、静水圧増分３Ｋε_ｍが計算される。その後、Ｓ５０５において、前回の計算ステップにおける各成分ごとの応力σ_ｉｊ ^０と、静水圧増分３Ｋε_ｍとに基づき、前述のようにして、更新された静水圧σ_ｍ ^ｔ＋ΔＴが計算される。

続いて、Ｓ５０６において、前記計算された各成分ごとのひずみ増分ε_ｉｊと等方ひずみ増分ε_ｍとに基づき、各成分ごとに偏差ひずみ増分ｅ_ｉｊが計算される。続いて、Ｓ５０７において、平面応力状態で、かつ、弾性域内において、偏差応力増分２Ｇｅ_ｉｊが計算される。

その後、Ｓ５０８において、前回の計算ステップにおける各成分ごとの偏差応力Ｓ_ｉｊ ^０と、偏差応力増分２Ｇｅ_ｉｊとに基づき、前述のようにして、更新された試行応力Ｓ_ｉｊ ^ｔ＋ΔＴが計算される。続いて、Ｓ５０９において、試行応力Ｓ_ｉｊ ^ｔ＋ΔＴの計算値に基づき、前述のようにして、相当応力の試行解σ^{Ｔｒｉａｌ}が計算される。

その後、図１４のＳ５１０において、板材がそれの両面の少なくとも一方において降伏したか否かが判定される。今回は、降伏したと仮定すれば、判定がＹＥＳとなり、Ｓ５１１以下のステップに移行する。

Ｓ５１１においては、図１５の式（１６）を用いることにより、相当塑性ひずみ増分ε^Ｐの増分Δε^Ｐが計算される。その後、Ｓ５１２において、その計算された増分Δε^Ｐに対応する相当応力σが対応相当応力σとされるとともに、その対応相当応力σと試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}とに対応する最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}が、図１５の式（１４）を用いて計算される。本実施形態においては、その最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}が前記（１）項における「暫定値」の一例である。

なお付言すれば、前記フォン・ミーゼスの降伏条件を採用しない場合には、前記ラジアルリターンアルゴリズムを、Ｓ５１１およびＳ５１２の反復的実行により実現することが可能である。

例えば、Ｓ５１１の各回の実行時には、相当塑性ひずみ増分ε^Ｐが一定の微小増分Δε^Ｐで増加させられる。これに対し、Ｓ５１２の各回の実行時には、相当塑性ひずみ増分ε^Ｐの暫定値（Ｓ５１０の実行によって計算される）に対応する相当応力σが対応相当応力σとして計算される。さらに、その対応相当応力σと試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}とに対応する最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値が計算される。

このＳ５１２においては、さらにまた、対応相当応力σと最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とが互いに十分に一致するか否かが判定される。十分に一致する場合には、相当塑性ひずみ増分ε^Ｐの暫定値と最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とがそれぞれ最終値とされる。最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の最終値は、更新された偏差応力成分を意味する。これに対し、対応相当応力σと最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}の暫定値とが互いに十分に一致しない場合には、Ｓ５１１に戻り、相当塑性ひずみ増分ε^Ｐの暫定値が一定の微小増分Δε^Ｐで増加させられる。

いずれの場合にも、前記ラジアルリターンアルゴリズムによって相当塑性ひずみ増分ε^Ｐの最終値と、更新された偏差応力成分Ｓ_ｉｊ ^{Ｆｉｎａｌ}とが取得されたならば、その後、Ｓ５１３において、更新された偏差応力成分Ｓ_ｉｊ ^{Ｆｉｎａｌ}と、前記更新された静水圧成分σ_ｍ ^ｔ＋Δｔとに基づき、図１５の式（１５）を用いることにより、各成分ごとに更新された応力σ_ｉｊ ^ｔ＋Δｔが計算される。

以上、板材がそれの両面の少なくとも一方において降伏した場合を説明したが、いずれの面においても降伏しない場合には、Ｓ５１０の判定がＮＯとなり、Ｓ５１４に移行する。このＳ５１４においては、前回の計算ステップまでに更新された試行応力Ｓ_ｉｊ ^{Ｔｒｉａｌ}と、前記更新された静水圧成分σ_ｍ ^ｔ＋Δｔとに基づき、図１５の式（１７）を用いることにより、各成分ごとに更新された応力σ_ｉｊ ^ｔ＋Δｔが計算される。

板材がそれの両面の少なくとも一方において降伏した場合には、Ｓ５１３の実行後、Ｓ５１５において、最終解Ｓ_ｉｊ ^{Ｆｉｎａｌ}（平面応力状態における面内応力σｘに相当する。）または前回の計算ステップにおける応力と、例えば図２８に示す式（４１）ないし図３０に示す式（５５）と、例えば図３３に示す式（７０）と、例えば図３３に示す式（７３）または（７５）とを用いることにより、前述のようにして、板厚方向応力σｚ（＝Ｆ（σｚｂ，σｚｃ)）を近似的に定義する応力近似関数が特定される。その特定される応力近似関数は、具体的には、曲げ変形に伴う板厚方向応力σｚｂの板厚方向における分布を定義する応力近似関数と、接触圧ｐによる板厚方向応力σｚｃの板厚方向における分布を定義する応力近似関数とを含んでいる。

その後、Ｓ５１６において、そのようにして特定された応力近似関数を用いることにより、前記ラジアルリターンアルゴリズムのもと、今回の板厚方向応力σｚを満たす曲げ方向応力σｘ’と断面方向応力σｙ’とが計算される。それら曲げ方向応力σｘ’と断面方向応力σｙ’とは、板厚方向における各位置ごとに、すなわち、前述の各積分点ごとに計算される。本実施形態においては、それら計算された曲げ方向応力σｘ’と断面方向応力σｙ’とがそれぞれ前記（１）項における「仮面内応力」の一例である。

続いて、Ｓ５１７において、上記特定された応力近似関数と、図３６において式（９１）で表す計算式とを用いることにより、前述のようにして、板厚方向応力σｚを考慮することに伴う断面力増分Ｆｚが計算される。この断面力増分Ｆｚは、各シェル要素ごとに計算される。そのシェル要素が面内積分点を１つしか有しない場合には、この断面力増分Ｆｚは、各シェル要素ごとに１つずつ計算される。

その後、Ｓ５１８において、その計算された断面力増分Ｆｚと、図３６において式（９２）で表す計算式とを用いることにより、前述のようにして、配分量σｚｍが計算される。この配分量σｚｍも、各シェル要素ごとに計算される。

続いて、Ｓ５１９において、その計算された配分量σｚｍと、図３６において式（９３）で表す計算式とを用いることにより、前述のようにして、板厚方向における各位置ごとに（例えば、各積分点ごとに）、曲げ方向応力σｘ’が補正されて曲げ方向応力の最終値σｘ’’が取得される。このＳ５１９においては、さらに、上記計算された配分量σｚｍと、図３６において式（９４）で表す計算式とを用いることにより、前述のようにして、板厚方向における各位置ごとに（例えば、各積分点ごとに）、断面方向応力σｙ’が補正されて断面方向応力の最終値σｙ’’が取得される。

すなわち、本実施形態においては、それら取得された最終値σｘ’’，σｙ’’がそれぞれ前記（１）項における「最終値」の一例を構成しているのである。また、図１２におけるＳ４０９において一連の解析が終了したと判定されたときにおけるそれら曲げ方向応力σｘと断面方向応力σｙとの最終値σｘ’’，σｙ’’がそれぞれ前記残留応力を意味している。その残留応力は、各シェル要素に関連付けてメモリ２２に記憶される。

その後、Ｓ５２０において、前回の計算ステップまでの各成分のひずみ増分と、今回の偏差ひずみ増分とに基づき、図１５の式（１９）を用いることにより、更新された塑性ひずみε^Ｐ _ｉｊ ^ｔ＋Δｔが計算される。この式（１９）における「Δλ」は、式（１８）および式（１９）から求められる。図１２におけるＳ４０９において一連の解析が終了したと判定されたときにおけるその塑性ひずみε^Ｐ _ｉｊ ^ｔ＋Δｔが前記残留ひずみを意味している。その残留ひずみは、各シェル要素に関連付けてメモリ２２に記憶される。

続いて、Ｓ５２１において、その更新された塑性ひずみε^Ｐ _ｉｊ ^ｔ＋Δｔに基づき、今回の注目シェル要素の板厚ｔ_０が更新される。この更新値は、応力近似関数の特定、ひずみ近似関数の特定、スプリングバックの解析等に影響を及ぼす。

いずれの場合にも、Ｓ５２２において、Ｓ５０１ないしＳ５２１の実行が板材のすべてのシェル要素について行われたか否かが判定される。Ｓ５０１ないしＳ５２１の実行がすべてのシェル要素について行われてはいないと仮定すれば、判定がＮＯとなり、Ｓ５２３において、次のシェル要素が新たな注目シェル要素とされた後、Ｓ５０１に戻る。これに対し、Ｓ５０１ないしＳ５２１の実行がすべてのシェル要素について行われたとすれば、Ｓ５２２の判定がＹＥＳとなり、Ｓ５２４において、板材を近似するシェル要素モデルに対して有限要素法が実行される。

具体的には、各シェル要素ごとの要素剛性マトリクスによって全体剛性マトリクスが構築され、各節点における外力と内力とのつり合いを考慮しつつ、連立の剛性方程式が解かれる。

その後、Ｓ５２５において、その剛性方程式の解が収束したか否かが判定される。収束しない場合には、判定がＮＯとなり、図１３のＳ５０１に戻り、新たな全体剛性マトリクスを用いてＳ５０１ないしＳ５２４が実行される。これに対し、上記剛性方程式の解が収束した場合には、Ｓ５２５の判定がＹＥＳとなり、以上で、この要素応力計算ルーチンの一回の実行が終了する。

その後、図１２のＳ４０４において、前記境界条件のもと、板材の変形とポンチ１１０の形状とを考慮することにより、ポンチ１１０から板材に作用する接触力が計算される。

続いて、Ｓ４０５において、ポンチ１１０の加工速度に関する境界条件を考慮することにより、各シェル要素の各節点の加速度が更新される。その後、Ｓ４０６において、ポンチ１１０と板材との接触に関する解析条件を考慮することにより、各シェル要素の各節点の加速度が更新される。続いて、Ｓ４０７において、それら節点加速度に対して時間積分を行うことにより、各節点の速度が更新され、さらに、その計算された速度に対して時間積分を行うことにより、各節点の変位が更新される。

その後、Ｓ４０８において、解析結果が出力されるとともに、必要なポスト処理が行われる。続いて、Ｓ４０９において、今回の解析が終了したか否かが判定される。Ｓ４０２ないしＳ４１０から成るループの一回の実行は、増分計算法における１回の時間増分（１回の計算ステップ）に対応しており、現象的には、ポンチ１１０がダイス１１２内に微小距離進入する動作に対応する。

したがって、Ｓ４０９においては、上記ループの初回の実行開始時刻からの経過時間が、そのループの所定実行回数に対応する長さに達したか否かが判定されることにより、今回の解析が終了したか否かが判定される。現在時刻ｔにおいては、未だ解析が終了してはいないと仮定すれば、判定がＮＯとなり、Ｓ４１０において、更新された現在時刻ｔ＋Δｔ（Δｔ：時間増分）が計算される。その後、Ｓ４０２に戻る。

Ｓ４０２ないしＳ４１０から成るループが何回か繰り返された結果、今回の解析が終了すれば、Ｓ４０９の判定がＹＥＳとなり、この成形シミュレーション・プログラムの一回の実行が終了する。

次に、前記材料特性値決定プログラムの内容を詳述する。

図３９には、この材料特性値決定プログラムがフローチャートで示されている。まず、Ｓ２０１において、材料が仮想的に分割された複数のシェル要素に順に付与された番号ｊが１に設定される。次に、Ｓ２０２において、メモリ２２から、今回のシェル要素ｊに対応する初期ひずみε^Ｐ _０（ｊ）が読み出される。

その後、Ｓ２０３において、その今回の初期ひずみε^Ｐ _０（ｊ）が、前記複数の予ひずみ区分のいずれに属するかが判定され、さらに、その属すると判定された予ひずみ区分に対応する材料特性値と、今回のシェル要素の番号ｊとを互いに関連付けるデータがメモリ２２に格納される。したがって、前記スプリングバック解析プログラムは、そのデータを参照することにより、各シェル要素ごとに選択された材料特性値を用いることにより、材料のスプリングバックを解析することができる。

続いて、Ｓ２０４において、今回の初期ひずみε^Ｐ _０（ｊ）に対して必要な修正が行われる。具体的には、上記Ｓ２０３において選択された材料特性値に厳密に対応する予ひずみε^Ｐ _ｇが厳密に今回の初期ひずみε^Ｐ _０（ｊ）に一致するとは限らないという事情を考慮し、一致しない場合には、一致するように今回の初期ひずみε^Ｐ _０（ｊ）が修正される。さらに、今回のシェル要素の番号ｊと、修正された初期ひずみε^Ｐ _０（ｊ）とを互いに関連付けるデータがメモリ２２に格納される。

したがって、前記スプリングバック解析プログラムは、そのデータを参照することにより、各シェル要素ごとに修正された材料特性値を用いることにより、材料のスプリングバックを解析することができる。

その後、Ｓ２０５において、番号ｊが最大値ｊＭＡＸ以上であるか否か、すなわち、材料の複数のシェル要素についてすべて材料特性値が選択されたか否かが判定される。今回は、番号ｊが最大値ｊＭＡＸ以上ではないと仮定すれば、判定がＮＯとなり、Ｓ２０６において、番号ｊが１増加させられ、その後、Ｓ２０２に移行する。以下、Ｓ２０２ないしＳ２０６が、新たなシェル要素について実行される。

Ｓ２０２ないしＳ２０６の実行が何回か繰り返された結果、Ｓ２０５の判定がＹＥＳとなれば、この材料特性値決定プログラムの一回の実行が終了する。

次に、前記バックストレス計算プログラムの内容を詳述する。

図４０には、このバックストレス計算プログラムがフローチャートで示されている。まず、Ｓ３０１において、材料の前記複数のシェル要素に順に付与された番号ｋが１に設定される。次に、Ｓ３０２において、メモリ２２から、今回のシェル要素ｋに対応する初期応力−ひずみ関係が読み出される。この初期応力−ひずみ関係は、材料の今回のシェル要素ｋにプレス成形直後に残存する応力とひずみとの関係を意味しており、前記成形シミュレーション・プログラムによりシミュレートされたものである。

なお、材料を仮想的に複数のシェル要素に分割する仕方は、バックストレスの計算と前記材料特性値の決定とで異ならないが、説明の便宜上、材料特性値の決定においては各シェル要素の番号がｊで表されているに対して、バックストレスαの計算においてはｋで表されている。

その後、Ｓ３０３において、その読み出された初期応力−ひずみ関係に基づき、今回のシェル要素ｋのバックストレスαが計算される。さらに、その計算されたバックストレスαと今回のシェル要素の番号ｋとを互いに関連付けるデータがメモリ２２に格納される。したがって、前記スプリングバック解析プログラムは、そのデータを参照することにより、前記材料特性値のみならず、各シェル要素ごとに計算されたバックストレスαをも用いることにより、材料のスプリングバックを解析することができる。

ここで、バックストレスαの計算手法を説明する。

プレス成形シミュレーションにおいては、プレス成形されるべき材料が仮想的に複数のシェル要素に分割されるとともに、各シェル要素について平面応力状態が仮定される。したがって、各シェル要素においては、図４１において式（１０１）で示すように、板厚方向のすべての応力成分が無視される。この式（１０１）において「σ」は垂直応力を意味し、「τ」はせん断応力を意味している。

このように仮定された各シェル要素においては、偏差応力ベクトル（σ’xx，σ’yy，σ’xy）の方向と、求めるべきバックストレスαの偏差成分（α’xx，α’yy，α’xy）の方向とが同じであると仮定すると、式（１０２）が成立する。この式（１２）において「Ａ」は定数を意味する。

一方、偏差応力ベクトルσ’ijについては、式（１０３）が成立する。この式（１０３）において「σij」は応力（垂直応力σのみならずせん断応力τも表す）、「δij」はクロネッカーのデルタ、「σm 」は静水圧下における応力σijをそれぞれ意味する。応力σm は、式（１０４）により求めることができる。

さらに、バックストレスαに関しては、式（１０５）が成立する。したがって、上記定数Ａが、図４２の式（１０６）により求められる。この式（１０６）において「ε」は、垂直ひずみεとせん断ひずみγとを選択的に表す。具体的には、式（１０６）の右辺における「σ」が垂直応力σを表す場合には垂直ひずみεを表し、これに対し、せん断応力τを表す場合にはせん断ひずみγを表す。

一方、上記式（１０２）により、図４２の式（１０７）で表される関係が成立する。さらに、バックストレスαの偏差成分α’については、偏差応力ベクトルσ’と同様に、式（１０８）が成立する。この式（１０８）において「αij」はバックストレス、「δij」はクロネッカーのデルタ、「αm 」は静水圧下におけるバックストレスαijをそれぞれ意味する。応力αm は、式（１０９）により求めることができる。

したがって、それら式（１０６）ないし式（１０９）を用いることにより、バックストレスαが計算される。図４３には、求めるべきバックストレスαの相当値が応力σの相当値との関係においてグラフで示されている。

以上のようにしてバックストレスαが計算された後、Ｓ３０４において、番号ｋが最大値ｋＭＡＸ以上であるか否か、すなわち、材料の複数の要素についてすべてバックストレスαが計算されたか否かが判定される。今回は、番号ｋが最大値ｋＭＡＸ以上ではないと仮定すれば、判定がＮＯとなり、Ｓ３０５において、番号ｋが１増加させられ、その後、Ｓ３０２に移行する。以下、Ｓ３０２ないしＳ３０５が、新たな要素について実行される。

Ｓ３０２ないしＳ３０５の実行が何回か繰り返された結果、Ｓ３０４の判定がＹＥＳとなれば、このバックストレス計算プログラムの一回の実行が終了する。

本発明者らは、本実施形態に従うスプリングバック解析方法の効果を検証するため、ワークＷとしてのある板材に対し、図４４に示すプレス成形機を用いてハット曲げ加工を実験的に行った。この実験により、後述の実験値が取得された。

図４４に示すように、そのプレス成形機は、図３に示すプレス成形機１００に対してブランクホルダ１２６が追加されたものに相当している。ブランクホルダ１２６は、板材のうちの縁部であって、ダイス１１２の上面によって支持されている部分に、プレス成形中、しわが発生することを防止するために、ダイス１１２の上面に対向する位置に配置されている。このブランクホルダ１２６は、プレス成形中、ダイス１１２と共同して、板材のうちの縁部を挟んで保持する。

図４４には、それらダイス１１２、ワークＷおよびブランクホルダ１２６の他にポンチ１１０も示されている。図４４には、そのポンチ１１０が、初期位置にある状態で実線で示される一方、フルストローク位置にある状態で破線で示されている。図４４には、さらに、ダイス１１２の肩部の丸味半径であるダイス肩Ｒが「Ｒｄ」で示され、また、ポンチ１１０の肩部の丸味半径であるポンチ肩Ｒが「Ｒｐ」で示されている。

図４４に示すプレス成形機を用いることにより、ハット曲げ加工が次の条件で実施された。

（１）板材に関する条件
・板厚ｔ_０：２．３［ｍｍ］
・曲げ半径Ｒ：３［ｍｍ］

（２）プレス成形機に関する条件
・ダイス肩Ｒｄ：３［ｍｍ］
・ポンチ肩Ｒｐ：５［ｍｍ］
・肩クリアランス：０．２３［ｍｍ］

本発明者らは、そのハット曲げ加工後に板材に発生するスプリングバックを評価するファクタとして、その板材の縦壁部の開き量と曲率とを選択した。

図４５に示すように、開き量Ｂは、板材の縦壁部が、それの幅方向（Ｘ方向）に、スプリングバックが発生しない理想形状（ダイス１１２の内面形状）から反った反り量として定義されている。

具体的には、この開き量Ｂを計算するために、板材の曲げ外面上に、底面からその板材の軸方向（Ｚ方向）に１５［ｍｍ］隔たった注目点ＰＸと、同様にして８０［ｍｍ］隔たった注目点ＰＣとが設定される。さらに、注目点ＰＸとＰＣとをつなぐ直線（図４５において破線で示す。）に沿って、注目点ＰＸから６５［ｍｍ］隔たった注目点ＰＹが設定される。

この注目点ＰＹが、プレス中心線（図４５において一点鎖線で示す。）から幅方向（Ｘ方向）に隔たった距離ｄが測定される。その測定された距離ｄから、ダイス１１２の内面の、プレス中心線からの幅方向隔たりＡ（例えば、５２．５３［ｍｍ］）を差し引くことにより、開き量Ｂが計算される。

これに対し、曲率は、図４５に示すように、板材の曲げ外面をその板材の中心線を含む平面で切断した場合にその板材のうちの縦壁部によって形成される曲線の曲率を意味するように定義されている。

具体的には、この曲率を計算するために、板材のうちの縦壁部の外面上に、底面から軸方向（Ｚ方向）に２０［ｍｍ］、５０［ｍｍ］および５０［ｍｍ］それぞれ隔たった３つの注目点ＰＡ、ＰＢおよびＰＣ点が設定される。それら注目点ＰＡ、ＰＢおよびＰＣを通過する単一の真円が想定され、その真円の半径ρの逆数として、曲率が計算される。

本発明者らは、上述の実験品と等価な板材に対して同じ条件でハット曲げ加工を行った場合にその板材に発生することが予想されるスプリングバックを３種類の成形シミュレーションによって解析した。図４６には、それら３種類の成形シミュレーションによる３種類の解析結果が、実験値と共に、開き量という観点からグラフで表されている。また、図４７には、それら３種類の成形シミュレーションによる３種類の解析結果が、実験値と共に、曲率（１／曲率半径ρ）という観点からグラフで表されている。

それら３種類の成形シミュレーションは、本実施形態に従うシミュレーションであって、板材の中立面を移動させることなく、板厚方向応力σｚを考慮することに伴う曲げ方向および断面方向における力の不つり合いを是正するものを含んでいる。図４６および図４７には、この種類の成形シミュレーションの結果が「軸力補正」というラベルを付して示されている。

それら３種類の成形シミュレーションは、さらに、板厚方向応力σｚを考慮することに伴う曲げ方向および断面方向における力の不つり合いを板材の中立面を移動させることによって是正するものも含んでいる。図４６および図４７には、この種類の成形シミュレーションの結果が「中立面移動」というラベルを付して示されている。

それら３種類の成形シミュレーションは、さらに、板厚方向応力σｚを考慮することに伴う曲げ方向および断面方向における力の不つり合いを是正するしないものも含んでいる。図４６および図４７には、この種類の成形シミュレーションの結果が「軸力補正なし」というラベルを付して示されている。

図４６に示すように、本実施形態に従う成形シミュレーションによれば、開き量の実験値が十分に高い精度で再現されるとともに、その再現精度は、中立面移動を伴う成形シミュレーションと同程度である。

また、図４７に示すように、本実施形態に従う成形シミュレーションによれば、曲率の実験値が十分に高い精度で再現されるとともに、その再現精度は、中立面移動を伴う成形シミュレーションより高い。

以上の説明から明らかなように、本実施形態においては、図１４におけるＳ５１２が前記（１）項における「暫定値算出工程」の一例を構成し、Ｓ５１５が同項における「応力近似関数特定工程」の一例を構成し、Ｓ５１６が同項における「仮面内応力算出工程」の一例を構成し、Ｓ５１７が同項における「断面力増分算出工程」の一例を構成し、Ｓ５１８が同項における「配分量決定工程」の一例を構成し、Ｓ５１９が同項における「最終値決定工程」の一例を構成しているのである。

さらに、本実施形態においては、図１４におけるＳ５１９が前記（２）項における「最終値算出工程」の一例を構成し、Ｓ５１８が前記（３）項における「配分量決定工程」の一例および前記（４）項における「配分量決定工程」の一例を構成し、Ｓ５１７が前記（５）項における「断面力増分算出工程」の一例を構成しているのである。

さらに、本実施形態においては、図１２に示す成形シミュレーション・プログラムが前記（８）項に係る「プログラム」の一例を構成し、その成形シミュレーション・プログラムを記憶しているメモリ２２、ＣＤ−ＲＯＭ３０またはＭＯ３２がそれぞれ、前記（９）項に係る「記録媒体」の一例を構成しているのである。

以上、本発明の実施の形態のうちの一つを図面に基づいて詳細に説明したが、これは例示であり、前記［発明の開示］の欄に記載の態様を始めとして、当業者の知識に基づいて種々の変形、改良を施した他の形態で本発明を実施することが可能である。

本発明の一実施形態に従う応力解析方法を含むスプリングバック解析方法を実施するのに好適なスプリングバック解析装置を示す系統図である。 上記スプリングバック解析方法を示す工程図である。 板材をＵ字状に曲げることが可能なプレス成形機の一例を示す正面断面図である。 上記スプリングバック解析方法において応力−ひずみ関係の実験値が予ひずみごとに取得される様子を説明するためのグラフである。 図２におけるステップＳ４のスプリングバック解析に使用される材料硬化モデルおよびそのスプリングバック解析に必要な材料特性値を説明するためのグラフである。 上記スプリングバック解析に必要な材料特性値を説明するための別のグラフである。 上記スプリングバック解析に使用される材料硬化モデルを等方硬化モデルおよび実際の応力−ひずみ関係と対比しつつ説明するためのグラフである。 図２におけるステップＳ２の材料特性値取得を行うためにコンピュータにより実行される材料特性値取得プログラムの内容を概念的に表すフローチャートである。 図８の材料特性値取得プログラムにより材料特性値が取得される理論および原理を説明するためのいくつかの式を示す図である。 図８におけるＳ１０６の詳細を材料軟化係数算出ルーチンとして概念的に表すフローチャートである。 図８の材料特性値取得プログラムにより予ひずみと材料特性値との関係が離散的に取得される様子を表形式で示す図である。 図２におけるステップＳ４の成形シミュレーションを行うためにコンピュータにより実行される成形シミュレーション・プログラムの内容を概念的に表すフローチャートである。 図１２におけるＳ４０３の詳細を要素応力計算ルーチンとして概念的に表すフローチャートの一部である。 上記要素応力計算ルーチンを表すフローチャートの別の一部である。 図１３および図１４に示す要素応力計算ルーチンに利用される理論および原理を説明するためのいくつかの式を示す図である。 図１３および図１４に示す要素応力計算ルーチンに利用される理論および原理を説明するためのグラフである。 従来例のシェル要素に利用される変位−ひずみ関係式を示す図である。 図１７の変位−ひずみ関係式により表現されるひずみ分布を板材の曲げ変形と共に示す図である。 断面回転角θｙと板材の曲げ内面からの距離ｄｓとの関係が曲げの厳しさＬｂに応じて変化する様子を示すグラフである。 図２におけるステップＳ４の成形シミュレーションを実行するための成形シミュレーション・プログラムに組み込まれたひずみ近似関数特定ルーチンの内容を概念的に表すフローチャートである。 断面回転角θｙと板材の曲げ内面からの距離ｄｓとの関係についてソリッド要素モデルを用いて計算された値の精度を実測値と対比して検証するためのグラフである。 図２０のひずみ近似関数特定ルーチンの実行結果のいくつかの例を示すグラフである。 せん断ひずみγｘｚと板材の曲げ内面からの距離ｄｓとの関係についての３種類の計算値であってソリッド要素モデルを用いた計算値と上記実施形態のシェル要素モデルを用いた計算値と従来例のシェル要素モデルを用いた計算値とから成るものを示すグラフである。 上記実施形態のシェル要素を用いて計算されたせん断ひずみγｘｚの分布と従来例のシェル要素を用いて計算されたせん断ひずみγｘｚの分布とを板材の曲げ変形と共に示す図である。 上記実施形態のシェル要素モデルに利用される変位−ひずみ関係式を示す図である。 上記実施形態のシェル要素を用いて計算された垂直ひずみεｘおよびせん断ひずみγｘｚの板厚方向における分布を板材の曲げ変形と共に示す図である。 上記実施形態のシェル要素を用いて計算された垂直応力σｚの分布と従来例のシェル要素を用いて計算された垂直応力σｚの分布とを板材の曲げ変形と共に示す図である。 図１４のＳ５１６において応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式の一部を示す図である。 上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式の別の一部を示す図である。 上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。 上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。 上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。 上記応力近似関数を特定するために利用される理論および原理を説明するためのいくつかの式のさらに別の一部を示す図である。 垂直応力σｚを考慮しない垂直応力σｘ，σｚと垂直応力σｚを考慮した垂直応力σｘ’，σｙ’との関係を説明するためのいくつかの式を示す図である。 図１４におけるＳ５１７ないしＳ５１９において、板材の中立面を移動させることなく、板厚方向応力σｚを考慮することに伴う曲げ方向および断面方向における力の不つり合いを是正するために行われる軸力補正の原理を概念的に説明するための図である。 上記軸力補正の原理を説明するためのいくつかの式を示す図である。 上記軸力補正を説明するためにシェル要素を示す斜視図である。 上記軸力補正を説明するために、板厚方向応力σｚの近似値を示すとともに、曲げ方向応力σｘを、板厚方向応力σｘを考慮しない場合の計算値と、考慮した場合の計算値と、上記軸力補正を行った場合の計算値とを示すグラフである。 図２におけるステップＳ５の材料特性値決定を行うためにコンピュータにより実行される材料特性値決定プログラムの内容を概念的に表すフローチャートである。 図２におけるステップＳ６のバックストレス計算を行うためにコンピュータにより実行されるバックストレス計算プログラムを概念的に表すフローチャートである。 図４０のバックストレス計算プログラムによりバックストレスを計算するために利用される理論および原理を説明するためのいくつかの式の一部を示す図である。 上記バックストレスを計算するために利用される理論および原理を説明するためのいくつかの式の別の一部を示す図である。 図４０のバックストレス計算プログラムの実行内容を説明するためのグラフである。 上記実施形態の効果を説明するためにハット曲げ加工を実施するために使用されたプレス成形機の要部を示す側面断面図である。 上記ハット曲げ加工によって板材に発生するスプリングバックを評価するためのファクタであるその板材の縦壁部の開き量および曲率を説明するための側面断面図である。 上記実施形態によるスプリングバックの解析精度を、図４５における開き量に関して、３つの比較例と共に説明するためのグラフである。 上記実施形態によるスプリングバックの解析精度を、図４５における曲率に関して、３つの比較例と共に説明するためのグラフである。

符号の説明

１２ コンピュータ
１６ 外部記憶装置
３０ ＣＤ−ＲＯＭ
３２ ＭＯ
１００ プレス成形機
１１０ ポンチ
１１２ ダイス
１１４ パッド
１２６ ブランクホルダ

Claims (9)

1. プレス成形による板材の曲げ変形をコンピュータによって解析するために、前記板材がそれの面方向に並んだ複数のシェル要素の集合体として前記コンピュータ上でモデル化されたシェル要素モデルを用いることにより、前記各シェル要素の板厚方向における応力である板厚方向応力を考慮した上で各シェル要素の面方向における応力である面内応力を解析する応力解析方法であって、
   前記各シェル要素につき、前記板厚方向応力が０である平面応力状態を想定して前記面内応力の暫定値を算出する暫定値算出工程と、
   その算出された面内応力の暫定値に基づき、前記板厚方向応力の前記板厚方向における分布を近似的に定義する応力近似関数を特定する応力近似関数特定工程と、
   その特定された応力近似関数を用いることにより、前記板厚方向における各位置ごとに、前記板厚方向応力と共に力学的に成立する面内応力を、前記板厚方向応力を考慮した仮面内応力として算出する仮面内応力算出工程と、
   前記板厚方向応力の前記板厚方向における分布に基づき、前記面内応力の算出に際して前記板厚方向応力を考慮することに起因して前記面内応力の、前記各シェル要素の断面上における総和である断面力が増加する量を断面力増分として算出する断面力増分算出工程と、
   その算出された断面力増分を、前記板厚方向に並んだ複数の位置にそれぞれ配分する配分量を決定する配分量決定工程と、
   前記板厚方向における各位置ごとに、前記算出された仮面内応力を、前記決定された配分量で補正することにより、前記面内応力の最終値を算出する最終値算出工程と
   を含む応力解析方法。
2. 前記最終値算出工程は、前記板厚方向における各位置ごとに、前記算出された仮面内応力から、前記決定された配分量を引き算することにより、前記算出された仮面内応力を補正する請求項１に記載の応力解析方法。
3. 前記配分量決定工程は、前記算出された断面力増分が前記複数の位置にそれぞれ互いに実質的に均等に配分されるように前記配分量を決定する請求項１または２に記載の応力解析方法。
4. 前記配分量決定工程は、前記板厚方向における位置に依存しないように前記配分量を決定する請求項１ないし３のいずれかに記載の応力解析方法。
5. 前記断面力増分算出工程は、前記応力分布関数を前記板厚方向における位置に関して実質的に積分することにより、前記断面力増分を算出する請求項１ないし４のいずれかに記載の応力解析方法。
6. 前記プレス成形は、前記板材の曲げ半径Ｒ［ｍｍ］をその板材の板厚ｔ［ｍｍ］で割り算して得られる曲げ指標が３より小さい条件で実施される請求項１ないし５のいずれかに記載の応力解析方法。
7. 前記プレス成形は、厚板成形とヘミング加工とハット曲げ加工との少なくとも一つを含む請求項１ないし６のいずれかに記載の応力解析方法。
8. 請求項１ないし７のいずれかに記載の方法を実施するためにコンピュータにより実行されるプログラム。
9. 請求項８に記載のプログラムをコンピュータ読み取り可能に記録した記録媒体。

Images