Material Modeling of Hot-Rolled Steel Sheet Considering Differential Hardening and Hole Expansion Simulation

Abstract

The elastic-plastic deformation behavior of a 440 MPa hot-rolled steel sheet subjected to many linear stress paths is precisely measured using biaxial tensile tests with cruciform specimens (ISO 16842: 2014) and multiaxial tube expansion tests (Kuwabara and Sugawara, 2013) to determine appropriate material models for finite element analysis (FEA). It was found that the Yld2000-2d yield function (Barlat et al., 2003) correctly reproduces the contours of plastic work (CPW) and the directions of the plastic strain rates (DPSR). Differential hardening (DH) models are determined by changing the values of exponent and material parameters of the Yld2000-2d yield function as functions of reference plastic strain. Moreover, FEA of the hole expansion forming of the test material is performed. The DH model correctly predicts the minimum thickness position that matches the fracture position of the specimen in experiment.

1. 緒言

自動車の足回り部品は高い応力が繰り返し負荷されることから比較的厚肉な熱延鋼板が用いられているため，車体総重量の多くを占めていながら軽量化が遅れていた¹⁾。軽量化のためには高強度材による薄肉化が有効だが，高強度化に伴う成形性の悪化によりプレス成形時に伸びフランジ部での割れが頻発し問題となっている²⁾。そこで，成形シミュレーションにより成形不具合を事前に予測し，金型修正工数を低減する試みが推進されている。成形不具合の予測精度を向上させるためには，材料の塑性変形挙動を再現可能な高精度な材料モデルが必要である³⁾。

伸びフランジ成形の数値解析結果に及ぼす材料モデルの影響に関する研究はいくつか報告されているが^4–9)，いずれの研究においても解析に用いた材料モデルの妥当性が立証されていないので，計算結果の妥当性に疑義が残る。これまで筆者らの研究グループは，十字形試験片を用いた二軸引張試験法^10,11)および円管材に大ひずみ二軸応力を発生させる二軸バルジ試験法¹²⁾を考案し，板材に対して適切な材料モデルを決定する実験手法として活用してきた。そして，張出し成形¹³⁾，面ひずみ¹⁴⁾，伸びフランジ成形^15,16)，液圧バルジ成形¹⁷⁾の成形シミュレーションにおいて，適切な降伏関数の選択が成形不具合の予測精度向上に必要不可欠であることを立証してきた。最近，Suzuki¹⁸⁾は，固溶強化型440 MPa級熱延鋼板(板厚2.0 mm)と析出強化型590 MPa級熱延鋼板(板厚2.3 mm)について穴広げ成形実験とその有限要素解析を実施し，高次降伏関数の方がHillの2次降伏関数よりも板厚分布の予測精度に優れることを示している。この結果は筆者らの研究グループにおいて得られた上記の知見と整合する。

成形シミュレーションでは，降伏曲面が相似形状を保ちつつ拡張すると仮定して，材料の加工硬化を定式化する材料モデル(等方硬化(Isotropic hardening:IH)モデル;以下IHモデル)が広く用いられている。しかし，実際の材料は加工硬化の進展に伴い降伏曲面の形状が変化することが知られている。これは異方硬化(Differential hardening:DH)^19,20)と呼ばれ，異方硬化を再現する材料モデル(以下DHモデル)²¹⁾を用いることで，成形シミュレーションの解析精度向上が期待できる。筆者らはこれまで，アルミニウム合金板^22–24)および鋼板²⁵⁾に対して，板材が破断するまでの大ひずみ範囲に対して二軸応力試験を実施し，Yld2000-2d降伏関数²⁶⁾を用いたDHモデルを作成して成形シミュレーションに適用することで，解析精度が向上することを示した。また，Tsutamoriら²⁷⁾は3次スプライン降伏関数を用いたDHモデルを作成し，穴広げ成形解析の精度が向上することを示している。

そこで本研究では，熱延鋼板を対象として，代表的な降伏関数に基づくIHモデルに加えて，パラメータの決定方法の異なる2種類のDHモデルを構築し，材料モデルが解析精度に及ぼす影響をより詳細に検証することを目的とする。研究の手法としては，供試材の単軸および二軸引張試験データに基づいて，IHおよびDHモデルを構築した。並行して，伸びフランジ成形を模擬する実験として穴広げ成形実験を行い，成形後の試験片の板厚分布を詳細に測定した。さらに，作成した材料モデルを用いて穴広げ成形解析を行い，板厚分布を計算し，実験値と比較した。その結果，DHモデルの適用は穴広げ成形における破断予測の精度向上に有効であることを明らかにしたので報告する。

2. 実験方法

2・1　供試材

供試材として公称板厚1.6 mmの440 MPa級熱間圧延鋼板を用いた。供試材の機械的性質をTable 1に示す。以下，異方性の主軸として，圧延方向(RD)をx軸，板幅方向(TD)をy軸，板厚方向(ND)をz軸とする。またRDから45°方向をDDと記す。

Table 1. Mechanical properties of the test material.

Tensile direction/º	E/GPa	c*/MPa	n*	α*	εpTS**	r***
0	192	764	0.215	0.012	0.187	0.78
45	203	756	0.234	0.023	0.191	1.06
90	212	797	0.249	0.029	0.182	0.85

* Approximated using σ=c(α+ε^p)ⁿ for ε^p=0.002~ε^p_TS

** Logarithmic plastic strain giving the maximum tensile load

*** Measured at a uniaxial nominal strain of ε_N=0.1

2・2　単軸引張試験

JIS13B号試験片を用いてRDから15°毎に単軸引張試験を行ない，r値および単軸塑性流動応力を測定した。以下，RDからθ°方向の単軸引張試験から測定されるr値をr_θ，単軸塑性流動応力をσ_θと記す。

2・3　二軸応力試験

供試材の二軸応力状態における弾塑性変形挙動を精度よく再現できる降伏関数を決定するために，十字形試験片による二軸引張試験^10,11)および円管試験片による二軸バルジ試験¹²⁾を行った。各試験片の寸法をFig.1に示す。

Fig. 1.

Geometry of the specimens of (a) biaxial tensile test (ISO 16842)⁴⁾ and (b) multiaxial tube expansion test⁶⁾. ↔ is the rolling direction of the original sheet sample. Units are in mm.

十字形試験片の形状および試験方法はISO16842¹⁰⁾に準拠した。Fig.1に示す位置においてひずみゲージによりひずみを測定することにより，応力の測定誤差を2%未満に抑えられることが有限要素解析^28,29)により確認されている。

二軸バルジ試験では，供試材をロール曲げして板縁を突合せ溶接した円管試験片を製作した。RDを円周方向にとった試験片IおよびRDを管軸方向にとった試験片IIの二種類を製作した。二軸バルジ試験法の詳細は文献12)を参照されたい。本試験法で測定された応力-ひずみ曲線は，円管試験片製作時の曲げによる予ひずみの影響を含むため，十字形試験片により測定した真応力-対数塑性ひずみ曲線により補正した。補正方法は文献12)を参照されたい。

両試験共に，σ_x:σ_y=4:1，2:1，4:3，1:1，3:4，1:2，1:4の7通りの線形応力経路にて実験を行った。ここで，σ_xはRD真応力，σ_yはTD真応力である。相当塑性ひずみ速度はほぼ10^-4 s^-1一定となるように制御した。二軸バルジ試験では，円管試験片の円周方向が常に最大主応力方向となるように試験を行った。すなわちσ_x:σ_y=4:1，2:1，4:3では試験片Iを，σ_x:σ_y=3:4，1:2，1:4では試験片IIを用いた。

3. 二軸応力試験結果

単軸および二軸応力状態における供試材の加工硬化挙動を定量的に評価し，材料モデルの高精度化に資するために，等塑性仕事面^19,20)を測定した。まず，RD単軸引張試験により得られた真応力-対数塑性ひずみ曲線において，規定の対数塑性ひずみε₀^p(以下，基準塑性ひずみ)に達した瞬間における単軸引張真応力σ₀と，ε₀^pに達するまでに消費された単位体積あたりの塑性仕事W₀を求める。次に，二軸応力試験およびTD単軸引張試験から得られた真応力-対数塑性ひずみ曲線においてW₀と等量の塑性仕事が消費された時点の真応力(σ_x，σ_y)を求め，それらを主応力空間にプロットし，規定のε₀^pに対する等塑性仕事面を決定した。等塑性仕事面の測定結果をFig.2に示す。単軸引張試験を含め，最大でε₀^p=0.22までの等塑性仕事面を測定できた。

Fig. 2.

Measured stress points forming contours of plastic work for different levels of ε₀^p.

等塑性仕事面の形状変化(異方硬化)を定量的に評価するために，等塑性仕事面を構成する応力点の値を当該のε₀^pに対応するσ₀で除して得られた無次元等塑性仕事面を求めた。さらに，各応力経路における形状比l/l_0.02を求めた。ここでl_0.02はε₀^p=0.02に対応する(σ_x/σ₀，σ_y/σ₀)までの原点からの距離，lは規定のε₀^pに対応する(σ_x/σ₀，σ_y/σ₀)までの原点からの距離である。二つの等塑性仕事面が相似であれば，すべての応力経路で形状比は1となる。ε₀^pの増加に伴うl/l_0.02の変化をFig.3に示す。ε₀^p≤0.02の範囲では材料の降伏伸びの影響により等塑性仕事面の形状変化が不安定なため省略した。0.02≤ε₀^p≤0.05の範囲でl/l_0.02は最大で3%変化し，ε₀^p>0.05では概ね一定となる。以上より，本供試材の異方硬化挙動が定量的に評価できた。このような等塑性仕事面の形状変化はIHモデルでは表現することはできない。

Fig. 3.

Variation of the shape ratio, l/l_0.02, with the increase of ε₀^p.

4. 材料モデリング

3節で測定された実験値に基づいて，von Mises³⁰⁾，Hillの2次(以下Hill’48)³¹⁾およびYld2000-2d²⁶⁾降伏関数(以下Yld2000-2d)を同定し，IHモデルを作成した。さらにYld2000-2dを用いてDHモデルを作成した。本節では各材料モデルの作成方法について述べる。

4・1　IHモデル

4・1・1　von MisesおよびHill’48

von Misesの同定にはσ₀を，Hill’48の同定には，σ₀および公称ひずみ0.10におけるr値r₀，r₄₅，r₉₀(Table 1参照)を用いた。以下，本モデルをr-Hill’48と記す。

4・1・2　Yld2000-2d

Yld2000-2dは次数Mおよび8つの未知係数α_i(i=1~8)より決定される。IHモデルの作成では，遺伝的アルゴリズム(Genetic algorithm:GA法)を用いて次式の評価関数fを最小化するMおよびα_i(i=1~8)を決定した。

  

$$f={\displaystyle \sum _{i=1}^N{w}_{\sigma ,i}{\left(l_{\text{M},i}-l_{\text{C},i}\right)}^2}+{\displaystyle \sum _{i=1}^N{w}_{\beta ,i}{\left({\beta }_{\text{M},i}-{\beta }_{\text{C},i}\right)}^2}$$(1)

第一項目はYld2000-2dから計算された降伏曲面と等塑性仕事面(実験値)との偏差を，第二項目は塑性ひずみ速度D^pの方向の計算値と実験値との偏差の総和を表す。具体的には，l_M,iおよびl_C,iは，それぞれi番目の線形応力経路における原点から無次元化等塑性仕事面および降伏曲面までの距離を，β_M,iおよびβ_C,i(単位:°)は，それぞれi番目の線形応力経路におけるD^pの方向の実験値とYld2000-2dによる計算値を表す(Fig.4参照)。w_σ,i，w_β,iは各応力経路に対する重み係数であり，この係数を変更することで，再現性を重視する実験値を任意に選択することが可能である。本研究において用いたw_σ,i，w_β,iの値をTable 2に示す。N(=14)は応力経路の数であり，7方向(RDから15°毎)の単軸引張試験および7線形応力経路の二軸引張試験から得られた実験値を用いた。等塑性仕事面を構成する応力の実験値については，単軸および二軸引張試験すべてにおいて，ε₀^p=0.22における塑性流動応力を用いた。D^pの方向の実験値については，単軸引張試験に対しては公称ひずみ10%におけるr値を，二軸引張試験に対してはε₀^p=0.22における実験値を用いた。

Fig. 4.

Schematic illustrations for identifying the parameter of the Yld2000-2d yield function: (a) contour of plastic work and (b) direction of plastic strain rate. φ_i is the angle of the i^th linear stress path from the RD in the σ_x-σ_y stress space.

Table 2. Weighting factors used in Eq.(1).

Stress ratio σx: σy		wσ,i	wβ,i
1: 0		1000	1
4: 1		1	5
2: 1		1	0.1
4: 3		1	0.1
1: 1		10	0.1
3: 4		1	0.1
1: 2		1	5
1: 4		1	0.1
0: 1		100	1
Uniaxial	15º	1	0.1
	30º	1	0.1
	45º	10	3
	60º	1	2
	75º	1	2

4・1・3　実験値とIHモデルによる計算値の比較

IHモデルにより計算された降伏曲面と等塑性仕事面の比較をFig.5(a)に，D^pの方向の比較をFig.5(b)に，σ_θの分布をFig.5(c)に，r_θの分布をFig.5(d)に示す。Yld2000-2dによる計算値が全ての実験値を最も精度よく再現した。

Fig. 5.

Measured data of (a) the stress points forming contours of plastic work, (b) the directions of plastic strain rates, (c) uniaxial flow stresses, and (d) r-values, compared with those calculated using selected yield functions (IH model).

4・2　DHモデル

DHモデルの作成では，ひずみの進展に伴うYld2000-2dのパラメータMおよびα_i(i=1~8)の変化を，ε₀^pの関数として近似した。Mおよびα_i(i=1~8)の発展を関数近似する方法としては，以下の二種類の方法(方法Aおよび方法B)を用いた。方法Aは文献²²⁾で考案された方法である。すなわち，規定のε₀^pにおけるMおよびα_i(i=1~8)を同定し，各パラメータを個別にε₀^pの関数として近似する。しかし，この方法ではMおよびα_i(i=1~8)がε₀^p毎に離散的に最適化されるため，ε₀^pの増加に伴う連続的な発展を再現しつつ最適化することが難しい。また，重み係数を変更する際に，その都度各ε₀^p毎のパラメータを同定する必要があるため，モデル作成に手間がかかる。そこで本報では新たに方法Bを考案した。方法Bでは，規定のε₀^pすべてに対して近似式の最適化を同時に行うことにより，連続的なパラメータの推移の最適化が可能となる。

4・2・1　方法Aに基づくDHモデル(DH-A)の作成

方法Aの手順を以下に記す。

i)規定のε₀^pにおけるMおよびα_i(i=1~8)を決定する。その決定手順は以下の通りである。次数Mを0.01刻みに変化させ，そのMに対して8つの実験値(規定のε₀^pにおけるσ₀，σ₄₅，σ₉₀および等二軸引張における流動応力σ_bと塑性ひずみ速度比 r_b≡ ${\dot{\epsilon }}_y^{\text{p}}/{\dot{\epsilon }}_x^{\text{p}}$ ，Table 1に示すr₀，r₄₅，r₉₀)を用いてNewton-Raphson法によりα_i(i=1~8)を決定し，降伏曲面と実験値の誤差を式(1)の評価関数fで計算し，fを最小化するMを決定する。ここでの重み係数は一律w_σ,i=1，w_β,i=0.3とした。

ii)i)において得られたMの変化をε₀^pの関数として次式で近似する。

  

$$M\left({\epsilon }_0^{\text{p}}\right)=\left(A-B\right)/\left[1+\exp \left\{\left({\epsilon }_0^{\text{p}}-C\right)/D\right\}\right]+B$$(2)

iii)式(2)で近似されたMを用いて，各ε₀^pに対するα_i(i=1~8)をGA法により決定する(計算方法は4・1・2節と同一)。

iv)iii)において得られたα_i(i=1~8)の変化をε₀^pの関数として次式で近似する。

  

$${\alpha }_i\left({\epsilon }_0^{\text{p}}\right)=A-B\exp \left(-C{\epsilon }_0^{\text{p}}\right)-D/\left(0.001+{\epsilon }_0^{\text{p}}\right)$$(3)

近似式(2)，(3)における係数A，B，C，Dの値はTable A1を参照。

4・2・2　方法Bに基づくDHモデル(DH-B)の作成

方法Bの手順を以下に記す。

i)4・2・1節と同様の手順でMおよびα_i(i=1~8)をε₀^pの関数として次式で近似する。

  

$$M\left({\epsilon }_0^{\text{p}}\right),{\alpha }_i\left({\epsilon }_0^{\text{p}}\right)\left(i=1~8\right)=A{\epsilon }_0^{\text{p}}\left({\epsilon }_0^{\text{p}}-B\right)\exp \left(C{\epsilon }_0^{\text{p}}\right)+D$$(4)

ii)式(4)の係数A，B，C，Dを初期値として，次式の評価関数f_Bが最小となるようにGA法を用いて最適化した。

  

$$f_{\text{B}}={\displaystyle \sum _{j=1}^m\left\{{\displaystyle \sum _{i=1}^n{w}_{\sigma ,i}{\left(l_{\text{M},i}-l_{\text{C},i}\right)}^2}+{\displaystyle \sum _{i=1}^n{w}_{\beta ,i}{\left({\beta }_{\text{M},i}-{\beta }_{\text{C},i}\right)}^2}\right\}}$$(5)

係数A，B，C，Dの数値はTable A2を参照。

4・2・3　実験値とDHモデルによる計算値の比較

DHモデルにより計算された降伏曲面とε₀^p=0.005および0.22における等塑性仕事面との比較をFig.6(a)に，D^pの比較をFig.6(b)に，σ_θの比較をFig.6(c)に，r_θの比較をFig.6(d)に示す。参考としてYld2000-2dのIHモデルによる計算値も併記する。DH-AおよびDH-Bに付随する表記「ε₀^p=X」は，ε₀^p=Xに対する実験値および計算値を表す。IHモデルはε₀^p=0.22における実験値を用いているため，ε₀^p=0.005における実験値と乖離が見られる。DHモデルは，ひずみの増大に伴う等塑性仕事面とσ_θの変化を精度よく再現している。r_θについては，DH-A，DH-Bともに実験値から若干の乖離があるものの，DDのr値が極大となる傾向は定性的に再現できている。

Fig. 6.

Measured data of (a) the stress points forming contours of plastic work, (b) the directions of plastic strain rates, (c) uniaxial flow stresses, and (d) r-values, compared with those calculated using the Yld2000-2d yield function (DH model).

5. 穴広げ成形実験と有限要素法解析

4節で作成した材料モデルを穴広げ成形シミュレーションに適用し，板厚ひずみ分布の測定値と比較することにより，作成した材料モデルの有用性を検証する。

5・1　穴広げ成形実験

穴広げ成形実験用金型の断面図をFig.7(a)に示す。直径215 mmの円形素板を使用し，ワイヤ放電加工により素板中央に直径30 mmの穴を開けた。パンチ径は100 mm，パンチ肩およびダイ肩の丸み半径は15 mmである。パンチストローク速度は0.1 mm/sとした。素板には，円周方向に10°，半径方向に2 mmごとの格子をけがき，試験後の板厚測定における初期座標の特定に用いた。パンチと試験片の間には両面にワセリンを塗布したテフロンシートを挿入し，潤滑した。パンチストロークhは28 mmとした。

Fig. 7.

Hole expansion forming: (a) geometry of the experimental apparatus and (b) quarter-model for finite element analysis (FEA).

5・2　穴広げ成形の有限要素解析

有限要素解析には，静的陰解法ソフトウェアABAQUS/Standard 2017を用いた。解析モデルの概要をFig.7(b)に示す。金型は解析的剛体とし，素板には4節点低減積分シェル要素S4Rを用い，板厚方向の積分点数は7とした。素板は対称性を考慮して1/4モデルで解析を行った。素板の要素は円周方向に3°毎，半径方向は1 mm毎に分割した。素板の初期板厚は実験に用いた試験片の平均初期板厚である1.574 mmとし，初期穴径は30 mmとした。ビード頂点部からの材料流入はないものとして，半径97.5 mmの素板外縁上の節点変位を完全固定とした。摩擦係数は素板とパンチ界面を0.03，金型と素板界面を0.2とした。また，加工硬化式にはRDの単軸引張試験から同定されたSwiftの式(Table 1参照)を用いた。

5・3　計算値と実験値の比較

5・3・1　IHモデル

初期素板上において穴縁から半径方向へ2 mm離れた位置における対数塑性板厚ひずみε_z^pの円周方向分布について，IHモデルによる計算値と実験値の比較をFig.8に示す。板厚ひずみの実験値はDDで極大，RD，TDで極小，RDで最小となった。等方性を仮定しているvon Misesでは当然ながら均一な板厚減少を示している。r-Hill’48は極大値と極小値の位相が実験値と逆転しており，材料モデルとしての精度は低い。これに対してYld2000-2dは極値の位置は実験値と一致している。しかしTDで最小値を予測しており，RDで最小となる実験の傾向を再現できていない。

Fig. 8.

Thickness strains along the expanded hole edge at a punch stroke of 28 mm, compared with the FEA results (IH model).

RD，DD，TDの半径方向のε_z^p分布について，IHモデルによる計算値と実験値の比較をFig.9に示す。Yld2000-2dはDDにおいて実験値とは最大で|Δε_z^p|=0.05の差異があるが，RDとTDでは実験値との差異は|Δε_z^p|=0.02以内であり，実験値の傾向を全般的に精度よく再現できている。またRD，DDにおいて穴縁から4~6 mm離れた位置で極小値をとる傾向も再現できている。

Fig. 9.

Measured and calculated thickness strains along RD, DD, and TD at a punch stroke of 28 mm (IH model).

5・3・2　DHモデル

円周方向ε_z^p分布について，DHモデルによる計算値と実験値の比較をFig.10に示す。DHモデルはDH-A，-BともにRDで最小，TDで極小となりかつ実験値ともほぼ一致しており，IHモデルよりも解析精度が向上したと言える。一方，DH-A，-Bどちらのモデルもε_z^pの振幅を過小評価している。またDDにおいてはΔε_z^p=0.05~0.04程度板厚減少を過大に予測している。DH-AとDH-Bの比較では，DH-Bの方がわずかながら実験値に近い。

Fig. 10.

Thickness strains along the expanded hole edge at a punch stroke of 28 mm, compared with the FEA results (DH model).

半径方向のε_z^p分布について，DHモデルによる計算値と実験値の比較をFig.11に示す。RD，DD，TDいずれの方向においても，3つの材料モデル(IH，DH-A，DH-B)間の差異はΔε_z^p≈0.01であり，実験値のばらつきと同程度である。RDとTDでは，どの材料モデルもΔε_z^p≈0.005~0.03の誤差範囲で実験値の傾向を再現できている。RDではDH-Aが，TDではIHとDH-Bが実験値に最も近い。一方，DDではIHが実験値に最も近いものの，いずれの材料モデルも0.02≤Δε_z^p≤0.06の範囲で板厚減少を過大に予測しており，RD，TDよりも実験値からの乖離が大きい。以上のように，半径方向のε_z^p分布の予測精度については材料モデル間の優劣はつけ難い。

Fig. 11.

Measured and calculated thickness strains along RD, DD, and TD at a punch stroke of 28 mm (DH model).

5・4　考察

穴広げ成形シミュレーションの結果を統括する。IHモデル同士の比較では，Yld2000-2dがε_z^p分布の実験値を最も精度よく再現した。さらにDHを再現するモデルの導入により，穴縁におけるε_z^pの最小位置を正確に再現することが出来た。これは，DHモデルではε₀^pの増加に伴うσ_θの円周方向分布の発展が的確に再現されたためと考える。σ_θの実験値とDH-Bによる計算値との比較をFig.12に示す。0.005≤ε₀^p≤0.01の低ひずみ範囲ではRDからTDにかけてσ_θは単調に大きくなるが，0.02≤ε₀^pではσ_θはDDで最小となった。DH-Bは実験値の傾向を精度よく再現できている。穴縁は単軸引張応力状態であるので，σ_θの発展の再現精度に優れたDHモデルがIHモデルよりも穴縁近傍のε_z^pの発達の予測精度が向上したものと推測する。実際，穴広げ試験においてさらにパンチを上昇させたところ，Fig.13に示すように，h=32 mmにおいてRDの穴縁より内側の位置で素板が破断した。これより，RDにおいてε_z^pが最小値をとることを的確に再現したDHモデルは伸びフランジ成形の破断予測に有効であると考える。この結果は文献18)および24)で得られた知見とも整合する。

Fig. 12.

Variation of σ_θ/σ₀ with the increase of ε₀^p, compared with that calculated using the DH-B model.

Fig. 13.

Top view of a specimen with the fracture in the RD at the hole edge.

一方，DHモデルをもってしても，DD半径方向の板厚減少が0.02≤Δε_z^p≤0.06の範囲で過大に予測され，RD，TDよりも実験値との乖離が大きくなった(Fig.11(b))。これは作成した材料モデルが依然として完璧ではないことを意味する。本研究では，σ_x-σ_y平面応力空間における等塑性仕事面とD^pの方向，ならびに穴縁近傍の変形状態と応力状態を代表するr_θとσ_θの分布の再現精度を高めるようにYld2000-2dのパラメータを決定した。一方，材料モデルの高精度化のためにはσ_x-σ_y-σ_xy空間における3次元形状の降伏曲面およびD^pの発展をくまなく高精度に再現する必要があるが，上記の材料試験のみでは，降伏曲面の発展を十分に捕捉できていない可能性がある。その対策として，以下に述べるような材料試験を追加で実施することが有効と考える。本供試材の異方性は比較的小さいことから，議論の簡単化のために等方性材料を仮定する。このとき，穴広げ成形においてDDに沿う応力状態は，σ_x-σ_y-σ_xy空間ではσ_x=σ_y平面上に載る。これらの応力状態を再現するためには，十字形試験片の腕方向をDDに平行に採り，応力比を変化させて二軸引張試験を行えばよい。こうして得られたDDの塑性変形特性をも再現できる材料モデルを構築できれば，DDの板厚分布の計算値はより実験値に近づく可能性がある。σ_x-σ_y-σ_xy空間における降伏曲面の発展を測定する他の方法としては，穴縁からある程度離れた位置における変形モードの代表として，応力の主軸方向を幾通りかに変化させて平面ひずみ引張試験(主ひずみ増分比 ${\dot{\epsilon }}_1^{\text{p}}:{\dot{\epsilon }}_2^{\text{p}}$ =1:0)³²⁾ を実施する方法も考えられる。

6. 結論

440 MPa級熱間圧延鋼板の二軸応力試験を行い，異方硬化挙動を精密に測定し，IHモデルに加えて異方硬化挙動を再現する2種類のDHモデルを作成した。作成した3種類の材料モデルを用いて穴広げ成形シミュレーションを行い，計算値と実験値を比較することにより，材料モデルが成形シミュレーションの解析精度に及ぼす影響を評価した。その結果，以下の知見を得た。

(1)σ_x-σ_y平面応力空間(σ_z=0)の第一象限に対して，9つの線形応力経路を付与した二軸応力試験による等塑性仕事面とD^pの方向の測定結果，およびRDから15°毎の単軸引張試験結果から，本供試材の異方硬化挙動を定量的に評価した。

(2)Yld2000-2dの次数Mおよび未知係数α_i(i=1~8)の発展をε₀^pの関数として表現することにより，DHモデルを作成した。

(3)(2)で作成したDHモデルを穴広げ成形シミュレーションに導入することで，RDにおいてε_z^pが最小となる実験値の傾向を的確に再現できた。ゆえに，伸びフランジ成形における破断予測の高精度化には，異方硬化挙動を再現できる材料モデルの使用が有効である。

(4)材料モデルのさらなる高精度化のためには，σ_x-σ_y-σ_xy空間における3次元形状の降伏曲面およびD^pの発展をくまなく高精度に再現する必要がある。そのためには，以下のような材料試験を追加実施し，得られた塑性変形特性をも再現できる材料モデルを構築することが有効と考える:(i)十字形試験片の腕方向をDDに平行に採り，応力比を幾通りかに変化させて二軸引張試験を行う;(ii)応力の主軸方向を幾通りかに変化させて平面ひずみ引張試験(主ひずみ増分比 ${\dot{\epsilon }}_1^{\text{p}}:{\dot{\epsilon }}_2^{\text{p}}$ =1:0)を実施する。

附録

DH-AおよびDH-Bモデルの近似式のパラメータ

4節にて作成したDH-Aモデルに用いた近似式(2)，(3)のパラメータをTable A1に，DH-Bモデルに用いた近似式(4)のパラメータをTable A2に示す。なお，DH-Aのα₅およびα₇に関してはひずみの進展による変化が小さかったため平均値を取り一定とした。

Table A1. Parameters of DH-A model.

	A	B	C	D
M*	6.295	8.1383	0.0521	0.0284
α1**	0.9882	0.1343	55.849	–0.0009
α2**	0.4932	–0.4932	0.0019	–0.0010
α3**	0.5137	–0.5141	–0.0233	0.0008
α4**	–5.6693	–6.682	–0.001	0.0004
α5**	1.0065	0		0
α6**	–15.318	–16.269	–0.0011	–0.0007
α7**	1.0145	0		0
α8**	1.0450	0.2079	108.236	0.0000

* Approximated using M(ɛ^p₀)=(A–B)/[1+exp{(ɛ^p₀–C)/D}]+B (Eq.(2))

** Approximated using α_i(ɛ^p₀)=A–B exp(–Cɛ^p₀)–D/(0.001+ɛ^p₀) (Eq.(3))

　

Table A2. Parameters of DH-B model.

	A	B	C	D
M*	–7.4392	–21.886	–53.627	7.3412
α1*	–74.845	3.6662	–574.46	1.0005
α2*	–8.4133	–38.792	–445.17	0.9623
α3*	–6.6872	–25.681	–344.42	0.9982
α4*	21.854	3.4405	–254.58	1.0098
α5*	2.6582	–6.4287	–380.71	1.0163
α6*	–147.99	2.0058	–747.27	1.0301
α7*	17.786	2.5921	–258.14	1.0238
α8*	29.91	4.3286	–378.23	1.0292

*Approximated using M, α_i(i=1~8)=Aɛ^p₀(ɛ^p₀–B) exp(Cɛ^p₀)+D (Eq.(4))

謝辞

供試材をご提供頂いた乃万暢賢博士(ユニプレス株式会社)に深甚なる謝意を表します。

文献

• 1)  M.Takahashi: J. Jpn. Soc. Technol. Plast., 58(2017), 105 (in Japanese). https://doi.org/10.9773/sosei.58.105
• 2)  E. Iizuka, T. Hira and A. Yoshitake: J. Jpn. Soc. Technol. Plast., 46(2005), 625 (in Japanese).
• 3)  F. Barlat and T. Kuwabara: J. Jpn. Soc. Technol. Plast., 57(2016), 230. https://doi.org/10.9773/sosei.57.230
• 4)  H.Kurosaki and Y.Unno: Trans. Jpn. Soc. Mech. Eng. C, 51(1985), 409 (in Japanese). https://doi.org/10.1299/kikaic.51.409
• 5)  H.Kurosaki, M.Tokiwa and K.Murai: Trans. Jpn. Soc. Mech. Eng. C, 52(1986), 380 (in Japanese). https://doi.org/10.1299/kikaic.52.380
• 6)  M.Gotoh, C.-R.Lim and M.Misawa: Trans. Jpn. Soc. Mech. Eng. C, 59(1993), 2855 (in Japanese). https://doi.org/10.1299/kikaic.59.2855
• 7)  H.Takuda, K.Mori, M.Kaneshiro and N.Hatta: Tetsu-to-Hagané, 84(1998), 182 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.84.3_182
• 8)  H.Takuda, T.Ozawa, T.Hama, T.Yoshida and J.Nitta: J. Jpn. Soc. Technol. Plast., 49(2008), 886 (in Japanese). https://doi.org/10.9773/sosei.49.886
• 9)  M.J. Worswick and M.J. Finn: Int. J. Plast., 16(2000), 701. https://doi.org/10.1016/S0749-6419(99)00069-8
• 10)  ISO 16842: 2014, Metallic materials –Sheet and strip –Biaxial tensile testing method using a cruciform test piece.
• 11)  T. Kuwabara, S. Ikeda and T. Kuroda: J. Mater. Process. Technol., 80-81(1998), 517. https://doi.org/10.1016/S0924-0136(98)00155-1
• 12)  T. Kuwabara and F. Sugawara: Int. J. Plast., 45(2013), 103. https://doi.org/10.1016/j.ijplas.2012.12.003
• 13)  T.Kuwabara, M.Umemura, K.Yoshida, M.Kuroda, S.Hirano and Y.Kikuta: J. Jpn. Inst. Light Met., 56(2006), 323 (in Japanese). https://doi.org/10.2464/jilm.56.323
• 14)  T.Moriya, T.Kuwabara, S.Kimura and S.Takahashi: J. Jpn. Soc. Technol. Plast., 51(2010), 43 (in Japanese). https://doi.org/10.9773/sosei.51.43
• 15)  T. Kuwabara, K. Hashimoto, E. Iizuka and J.W. Yoon: J. Mater. Process. Technol., 211(2011), 475. https://doi.org/10.1016/j.jmatprotec.2010.10.025
• 16)  K.Hashimoto, T.Kuwabara, E.Iizuka and J.W.Yoon: Tetsu-to-Hagané, 96(2010), 557 (in Japanese). https://doi.org/10.2355/tetsutohagane.96.557
• 17)  D. Yanaga, T. Kuwabara, N. Uema and M. Asano: Int. J. Solids Struct., 49(2012), 3488. https://doi.org/10.1016/j.ijsolstr.2012.03.005
• 18)  T.Suzuki: Ph.D. thesis, Hiroshima University, (2018), (in Japanese), http://ir.lib.hiroshima-u.ac.jp/00046059, (accessed 2020-03-12).
• 19)  R. Hill and J.W. Hutchinson: J. Appl. Mech., 59(1992), S1. https://doi.org/10.1115/1.2899489
• 20)  R. Hill, R., S.S.Hecker and M.G.Stout: Int. J. Solids Struct., 31(1994), 2999. https://doi.org/10.1016/0020-7683(94)90065-5
• 21)  T.Kuwabara, H.Takizawa and T.Hakoyama: J. Jpn. Soc. Technol. Plast., 57(2016), 181 (in Japanese). https://doi.org/10.9773/sosei.57.181
• 22)  T.Yanaga, H.Takizawa and T.Kuwabara: J. Jpn. Soc. Technol. Plast., 55(2014), 55 (in Japanese). https://doi.org/10.9773/sosei.55.55
• 23)  J.Kawaguchi T.Kuwabara and T.Sakurai: J. Jpn. Inst. Light Met., 65(2015), 554 (in Japanese). https://doi.org/10.2464/jilm.65.554
• 24)  T. Kuwabara, T. Mori, M. Asano, T. Hakoyama and F. Barlat: Int. J. Plast., 93(2017), 164. https://doi.org/10.1016/j.ijplas.2016.10.002
• 25)  T.Hakoyama and T.Kuwabara: From Creep Damage Mechanics to Homogenization Methods, ed. by H.Altenbach, T.Matsuda and D.Okumura, Springer, Heidelberg, (2015), 67. https://doi.org/10.1007/978-3-319-19440-0_4
• 26)  F. Barlat, J.C. Brem, J.W. Yoon, K. Chung, R.E. Dick, D.J. Lege, F. Pourboghrat, S.-H. Choi and E. Chu: Int. J. Plast., 19(2003), 1297. https://doi.org/10.1016/S0749-6419(02)00019-0
• 27)  H.Tsutamori, E.Iizuka, T.Amaishi, K.Sato, Y.Ogiwara and J.Miyamoto: J. Jpn. Soc. Technol. Plast., 57(2016), 245 (in Japanese). https://doi.org/10.9773/sosei.57.245
• 28)  Y. Hanabusa, H. Takizawa and T. Kuwabara: J. Mater. Process. Technol., 213(2013), 961. https://doi.org/10.1016/j.jmatprotec.2012.12.007
• 29)  Y. Hanabusa, H. Takizawa and T. Kuwabara: Steel Res. Int., 81(2010), 1376.
• 30)  R.v.Mises: Göttingen Nachrichten, Math.-phys. Klasse, (1913), 582 (in German). http://eudml.org/doc/58894, (accessed 2020-01-11).
• 31)  R. Hill: Proc. Roy. Soc. Lond. A, 193(1948), 281. https://doi.org/10.1098/rspa.1948.0045
• 32)  T.Kuwabara and S.Ikeda: Tetsu-to-Hagané, 88(2002), 334 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.88.6_334