JP4532143B2 - 板材成形シミュレーション及びプレス成形方法 - Google Patents

Description

本発明は、シェル要素を用いた有限要素法に基づく板材成形シミュレーション及びプレス成形方法に関する。

近年、板材のプレス成形分野において、有限要素法に基づく板材成形シミュレーションが普及してきた。有限要素法とは、解析対象となる板材を幾つかの要素に分割して考え、要素毎に方程式を作り、それを元に解析対象全体の方程式（全体方程式）を組み立てて解く解法である。

現在、板材成形シミュレーションにおいて広く利用されているシェル要素としては、Mindlin-Reissnerの仮定を用いたものが知られている（非特許文献１等参照）。このシェル要素は、板材の厚さ方向の応力分布を常に一定とした平面応力状態を仮定している。

Ｏ．Ｃ．ツェンキビッツ他、マトリックス有限要素法、科学技術出版社、１９９８年１月

従来の板材成形シミュレーションでは、板材の厚さ方向の応力を常に「ゼロ」と仮定したシェル要素を用いているため、厚さ方向のひずみ及び応力を求めることはできない。従って、板厚方向の応力を無視した状態でのシミュレーションを行うため、応力精度に敏感なスプリングバック変形予測の精度を悪化させていた。

また、従来の板材成形シミュレーションでは、厚さ方向に節点及び要素が一層しか存在しないシェル要素を用いているため、シェル要素の上下面に異なる境界条件を与えることができない。従って、板材の下面を拘束して、板材の上面に荷重を加えるような曲げ変形をシミュレートすることができなかった。

ところで、板材成形シミュレーションに、シェル要素に代えて、ソリッド要素を用いれば、厚さ方向のひずみ及び応力を求めること、及びソリッド要素の上下面に異なる境界条件を与えることができるが、計算時間がシェル要素を用いた場合の数倍になってしまう。そのため、計算時間の制約から、板材成形シミュレーションにソリッド要素を用いることはほとんどない。

そこで、本発明の目的は、板材の厚さ方向のひずみ及び応力を精確に求めることができる、シェル要素を用いた有限要素法に基づく板材成形シミュレーション及びプレス成形方法を提供することにある。

上記目的を達成するために、本発明は、四角形のシェル要素を用いた有限要素法に基づいて板材の曲げ変形をシミュレートする板材成形シミュレーション方法において、コンピュータに上記板材の形状が入力されると、上記コンピュータが、上記板材をシェル要素の四隅の節点に上下面中央の節点を加えた六節点のシェル要素で定義すると共に、そのシェル要素の形状関数を以下の［数１］、［数２］、［数３］で定義し、

（ただし、［数１］〜［数３］において、uはx方向の変位、vはy方向の変位、wはz方向（板厚方向）の変位、θxはx軸回りの回転角度、θyはy軸回りの回転角度、ξはx方向に対応した局所座標系により表される座標、ηはy方向に対応した局所座標系により表される座標、ζはz方向に対応した局所座標系により表される座標を示す）
上記コンピュータが、上記［数１］、［数２］、［数３］で定義される形状関数を用いて各節点の変位を算出すると共に、上記コンピュータが、上記［数１］、［数２］、［数３］で定義される形状関数を微分して上下面中央の節点のひずみを求め、求めたひずみから応力−ひずみ関係を通じて上下面中央の節点の応力を求め、上記コンピュータが、求めたひずみ及び応力から上記板材を曲げ又は上記板材を板厚方向に圧縮した際に発生する板厚方向の垂直方向変位を計算し、上記コンピュータが、計算した板厚方向の垂直方向変位から上記板材の曲げ変形に基づくスプリングバック量を算出することを特徴とする板材成形シミュレーション方法である。

本発明によれば、板材の厚さ方向のひずみ及び応力を精確に求めることができるという優れた効果を奏する。

以下、本発明の好適な一実施形態を添付図面に基づいて詳述する。

まず、図１を用いて、本実施の形態において用いるシェル要素について説明する。

図１は、本実施の形態に用いるシェル要素の概略図である。図１において、図中の横方向をx方向、図中の斜め方向をy方向、図中の上方向をz方向（板材の厚さ方向）とそれぞれ定義する。

図１に示すように、本実施の形態に用いるシェル要素ｅは、四角形の平面要素である。このシェル要素ｅのx方向の長さ、y方向の長さ、及びz方向の長さは、それぞれL1、L2、及びtであるとする。

シェル要素ｅの各隅には、四つの第一の節点ｎ１〜ｎ４が配置されている。シェル要素ｅの上下面中央には、二つの第二の節点ｎ５、ｎ６が配置されている。第二の節点ｎ５は、解析対象となる板材の下面に相当する。第二の節点ｎ６は、解析対象となる板材の上面に相当する。

シェル要素ｅの形状関数を数１〜数３に示す。

数１〜数３において、u、v及びwは、それぞれx、y、及びz方向の変位である。θ_x及びθ_yは、それぞれx軸、及びy軸回りの回転角度である。また、これらの記号に付された下付数字は、各節点ｎ１〜ｎ６の番号と対応している（例えば、u₁は、第一の節点ｎ１のx方向の変位である）。ξ、η、及びζは、それぞれx、y、及びz方向に対応した、局所座標系により表される座標を示す。

数１〜数３に示すように、節点ｎ１〜ｎ４は、x、y、z方向の変位、及びx軸、y軸回りの回転角度が定義されている。これら節点ｎ１〜ｎ４は、背景技術の欄で説明した平面応力状態を仮定している。また、節点ｎ５、ｎ６は、それぞれ、z方向（板材の厚さ方向）の変位のみが定義された一自由度の節点である。

つまり、本実施の形態に用いるシェル要素ｅは、四角形のシェル要素ｅの四隅の第一の節点ｎ１〜ｎ４に、上下面中央の第二の節点ｎ５、ｎ６を加えた六節点のシェル要素ｅである。

ここで、シェル要素ｅに、z方向（板材の厚さ方向）の変位が定義された第二の節点ｎ５、ｎ６を、板材の上下面に相当する位置にそれぞれ配置している。そのため、シェル要素ｅの上下面に異なる境界条件を与えることができると共に、第二の節点ｎ５、ｎ６の変位に基づいて、厚さ方向のひずみ及び応力を求めることができる。

また、シェル要素ｅは、従来のシェル要素に、一節点につき一自由度を有する節点（ｎ５、ｎ６）を二個加えたものである。そのため、従来のシェル要素に比べて計算時間が若干増加する程度であると共に、ソリッド要素に比べて計算時間が短いことを維持できる。

次に、数４及び数５に、シェル要素ｅの応力を定義する。

数４において、σ_xx、σ_yy、及びσ_zzは、それぞれ、シェル要素ｅ内のx、y、及びz方向の垂直応力である。σ_xy、σ_yz、及びσ_xzは、それぞれ、xy平面、yz平面、及びxz平面のせん断応力である。

従来のシェル要素では、z方向の垂直応力σ_zzを、常にゼロと仮定していた。本実施の形態の一例として、z方向の垂直応力σ_zzを、数５のように定義した。数５において、Ｐは、各節点ｎ１〜ｎ６に付与する荷重を示す。

これら数１〜数５を、予めコンピュータの記録媒体に記憶させておく。また、有限要素法に基づく板材成形シミュレーションに用いる数１〜数５以外の関数（剛性関数等）は、従来と同様のものが使用可能であるので、説明を省略する。これらもコンピュータの記憶媒体に予め記憶させておく。

次に、図２を用いて、本実施の形態の板材成形シミュレーションの手順について説明する。

まず、ステップＳ１において、コンピュータに解析対象となる板材の形状データ（板材の大きさ）、板材の物性値（ヤング率、ポアソン比等）、及び解析モデルの定義に使用するシェル要素ｅの数（分割数）等を入力する。次に、ステップＳ２において、板材の形状データ等に基づいて、解析モデルをシェル要素ｅで定義する。つまり、シェル要素ｅの形状関数を元に全体方程式を構築する。

その後、ステップＳ３において、コンピュータに境界条件及び荷重条件等を入力する。次に、ステップＳ４において、全体方程式を解くことにより、各節点ｎ１〜ｎ６の変位を算出すると共に、ステップＳ５において、第二の節点ｎ５、ｎ６のひずみ及び応力を求め、その求めた第二の節点ｎ５、ｎ６のひずみ及び応力から、板材の曲げ変形に基づくスプリングバック量を算出する。そして、ステップＳ６において、ステップ５で求めたスプリングバック量を出力する。

ここで、本実施の形態では、ステップＳ４において算出した各節点ｎ１〜ｎ６の変位に基づいて、境界条件及び荷重条件等を変更して再度シミュレートすることができる。この場合、ステップＳ２に戻り、ステップＳ４で算出した変位に基づいて、変形後の全体方程式を構築する。

次に、図３を用いて、本実施の形態の板材成形シミュレーションの一例について説明する。

図３（ａ）は、解析モデルの一例を示す概略図である。図３（ｂ）は、図３（ａ）の解析モデルを下方へたわませた状態を示す概略図である。

まず、図３（ａ）に示すように、成形すべき板材の解析モデルを、任意数（ここでは、二個）のシェル要素ｅ１、ｅ２で定義する。これらは、図２に示すステップＳ１「形状データの入力」、及びステップＳ２「全体方程式の構築」に相当する。ここで、節点ｍ１、ｍ２が板材の一端を構成し、節点ｍ７、ｍ８が板材の他端を構成する。節点ｍ３、ｍ４は、板材の長手方向の中央部を構成する。この節点ｍ３、ｍ４は、シェル要素ｅ１とシェル要素ｅ２とで共有される。また、節点ｍ５、ｍ９は、板材の下面を構成する。節点ｍ６、ｍ１０は、板材の上面を構成する。

これら二つのシェル要素ｅ１、ｅ２は、図１に示すシェル要素ｅと同様のものである。シェル要素ｅ１については、節点ｍ１〜ｍ４がシェル要素ｅの第一の節点ｎ１〜ｎ４、節点ｍ５、ｍ６がシェル要素ｅの第二の節点ｎ５、ｎ６にそれぞれ相当する。シェル要素ｅ２については、節点ｍ３、ｍ４、ｍ７、ｍ８がシェル要素ｅの第一の節点ｎ１〜ｎ４、節点ｍ９、ｍ１０がシェル要素ｅの第二の節点ｎ５、ｎ６にそれぞれ相当する。

次に、節点ｍ１、ｍ２、及び節点ｍ７、ｍ８にx、z方向への変位を拘束する境界条件を与えた後、節点ｍ１、ｍ２、及び節点ｍ７、ｍ８に回転増分Ａを与え、解析モデルを下方へたわませる。これは、図２に示すステップＳ３「境界条件の入力」に相当する。回転増分Ａの付与は、節点ｍ１、ｍ２、及び節点ｍ７、ｍ８をy軸回りにθ度だけ回転させることによって行われる。

このときの各節点ｍ１〜ｍ１０の変位を算出し、図３（ｂ）に示す解析モデルを構築する。これは、図２に示すステップＳ４「変位の算出」、及びステップＳ２「全体方程式の構築」に相当する。

次に、図３（ｂ）に示す解析モデルの節点ｍ６、ｍ１０にz方向に垂直に荷重Ｐを付与する。図４に示すものは、板材の下面に相当する節点ｍ５、ｍ９にz方向への変位を拘束する境界条件を与えないものであり、図５に示すものは、板材の下面に相当する節点ｍ５、ｍ９にz方向への変位を拘束する境界条件を与えたものである。

このときの各節点ｍ１〜ｍ１０の変位を算出すると共に、第二の節点ｍ５、ｍ６、ｍ９、ｍ１０のひずみ及び応力を求める。これは、図２に示すステップＳ４「変位の算出」に相当する。

図６に、図４及び図５の解析モデルにおけるz方向（板材の厚さ方向）の応力分布を示す。図６において、横軸は、z方向の応力σ_zzを示し、縦軸は、z方向の座標ζを表す。ζ＝１のときの応力σ_zzは、板材の上面、即ち節点ｍ６（ｍ１０）の応力を示す。ζ＝−１のときの応力σ_zzは、板材の下面、即ち節点ｍ５（ｍ９）の応力を示す。これらを繋いで、板材の厚さ方向の応力分布とする。

実線ｒ１及び破線ｒ２は、それぞれ、図４の解析モデルの応力分布、及びその理論的な応力分布を示す。実線ｒ３及び破線ｒ４は、それぞれ、図５の解析モデルの応力分布、及びその理論的な応力分布を示す。

ここで、図６に示すように、算出した応力分布ｒ１、ｒ３と理論的な応力分布ｒ２、ｒ４とでは、それぞれ若干の誤差が生じている。これは、節点ｍ５、ｍ６、ｍ９、ｍ１０（シェル要素ｅの第二の節点ｎ５、ｎ６）がz方向（板材の厚さ方向）の変位のみが定義されているためである。シェル要素ｅの第二の節点ｎ５、ｎ６に第一の節点ｎ１〜ｎ４と同様の自由度を付与すれば、算出した応力分布ｒ１、ｒ３と理論的な応力分布ｒ２、ｒ４とは完全に一致するが、当然、計算時間が数倍に増加する。

次に、算出した応力分布に基づいて、図５の解析モデルにおけるスプリングバック量の予測をする。

図７において、横軸は、荷重Ｐを示し、横軸は、解析モデルの長手方向端部のy軸回りの回転モーメント（曲げモーメント）Ｍ（図５参照）を示す。この曲げモーメントＭは、図６で示した応力分布より求めるものとする。

図７に示すように、荷重Ｐが増加すると、解析モデルの長手方向端部の曲げモーメントＭが低下している。この曲げモーメントＭは、板材成形におけるスプリングバックの原因となり、成形後の板材をスプリングバックさせる。本実施の形態においては、プレス成形時のスプリングバック量を、曲げモーメントＭの大きさから予測するものとする。

本実施の形態のプレス成形方法は、上述の板材成形シミュレーションにより求めたスプリングバック量を考慮して、実際に板材をプレス成形するものである。詳しくは、スプリングバック量を考慮して、成形圧力（荷重）及び金型の形状等を決定するものとする。

以上、本実施の形態の板材成形シミュレーションは、板材を、従来のシェル要素の上下面中央に節点ｎ５、ｎ６を加えた六節点のシェル要素で定義している。そのため、本実施の形態の板材成形シミュレーションによれば、ソリッド要素を用いることなく、板材の厚さ方向の応力分布を精確に求めることができる。これにより、板材の厚さ方向の応力分布から、板材の曲げ変形に基づくスプリングバック量を予測することができる。

また、本実施の形態のプレス成形方法は、上述の板材成形シミュレーションにより得られたスプリングバック量を考慮してプレス成形を行う。そのため、本実施の形態のプレス成形方法によれば、スプリングバック量を考慮した板材の変形を予測することができ、従来行っていた「決め押し」等の工程を行う必要がなくなる。これにより、従来に比べて加工時間（成形時間）を大幅に短縮することが可能となる。

本実施の形態に用いるシェル要素の概略図である。 板材成形シミュレーションのステップ図である。 （ａ）は、解析モデルの一例を示す概略図である。（ｂ）は、図３（ａ）の解析モデルを下方へたわませた状態を示す概略図である。 図３（ｂ）の解析モデルに荷重を付加した状態を示す概略図である。 図３（ｂ）の解析モデルに荷重を付加した状態を示す概略図である。 z方向の応力分布を示すグラフである。 荷重と解析モデルの長手方向端部の曲げモーメントとの関係を示すグラフである。

符号の説明

ｅ シェル要素
ｎ１、ｎ２、ｎ３、ｎ４ 第一の節点
ｎ５、ｎ６ 第二の節点
u、v、w 変位
θ 回転角度
σ 応力
Ｐ 荷重
Ｍ 曲げモーメント

Claims (1)

1. 四角形のシェル要素を用いた有限要素法に基づいて板材の曲げ変形をシミュレートする板材成形シミュレーション方法において、
   コンピュータに上記板材の形状が入力されると、
   上記コンピュータが、上記板材をシェル要素の四隅の節点に上下面中央の節点を加えた六節点のシェル要素で定義すると共に、そのシェル要素の形状関数を以下の［数１］、［数２］、［数３］で定義し、
   

   

   
   

   

   
   

   

   
   （ただし、［数１］〜［数３］において、uはx方向の変位、vはy方向の変位、wはz方向（板厚方向）の変位、θxはx軸回りの回転角度、θyはy軸回りの回転角度、ξはx方向に対応した局所座標系により表される座標、ηはy方向に対応した局所座標系により表される座標、ζはz方向に対応した局所座標系により表される座標を示す）
   上記コンピュータが、上記［数１］、［数２］、［数３］で定義される形状関数を用いて各節点の変位を算出すると共に、
   上記コンピュータが、上記［数１］、［数２］、［数３］で定義される形状関数を微分して上下面中央の節点のひずみを求め、求めたひずみから応力−ひずみ関係を通じて上下面中央の節点の応力を求め、
   上記コンピュータが、求めたひずみ及び応力から上記板材を曲げ又は上記板材を板厚方向に圧縮した際に発生する板厚方向の垂直方向変位を計算し、
   上記コンピュータが、計算した板厚方向の垂直方向変位から上記板材の曲げ変形に基づくスプリングバック量を算出することを特徴とする板材成形シミュレーション方法。

Images