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

Abstract

The elastoplastic 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 to determine appropriate material models for finite element analysis (FEA). It is found that the Yld2000-2d yield function correctly reproduces the contours of plastic work and the directions of the plastic strain rates. Differential hardening (DH) models are determined by varying the values of the exponent and material parameters for the Yld2000-2d yield function as functions of the reference plastic strain. Moreover, a finite element analysis of hole expansion in the test material is performed. The DH model correctly predicts the minimum thickness position, which matches the fracture position of the specimen in the experiment.

1. Introduction

Thick hot-rolled steel sheets are widely applied to the undercarriage parts of automobiles, which are subjected to repeated high stress. Therefore, weight reduction of these parts has been delayed even though they account for a large part of the total weight of the car body.¹⁾ Although thinning using a high-strength material effectively reduces weight, the deterioration of formability due to an increase in strength frequently leads to fracture in the stretch flange area during press forming.²⁾ Therefore, attempts have been made to predict fracture in advance via forming simulations and to reduce the effort required for establishing optimal die shape and forming conditions. To improve the accuracy of forming defect prediction, a highly accurate material model that can reproduce the plastic deformation behavior of a material during forming operations is required.³⁾

Several studies have investigated the effect of the material model on the numerical analysis results of stretch flange forming.^4,5,6,7,8,9) However, none proved the validity of the material model used in the analysis, and thus the validity of the calculation results is in doubt. Our research group devised a biaxial tensile test method^10,11) that uses a cruciform test piece and a multiaxial tube expansion test method¹²⁾ that generates a biaxial stress for a large strain range in a tubular or sheet material. These methods were used to experimentally determine a proper material model for tubular and sheet materials. It was demonstrated that the selection of an appropriate yield function improves the accuracy of defect prediction in forming simulations of stretching,¹³⁾ minute surface deflection,¹⁴⁾ stretch flange forming,^15,16) and hydraulic bulge forming.¹⁷⁾ Recently, Suzuki¹⁸⁾ conducted hole expansion forming experiments and a finite element analysis (FEA) of a solute-strengthened 440-MPa hot-rolled steel sheet (sheet thickness: 2.0 mm) and a precipitation-strengthened 590-MPa hot-rolled steel sheet (sheet thickness: 2.3 mm) and found that a higher-order yield function is superior to Hill’s quadratic yield function in predicting the thickness distribution. This result is consistent with the findings obtained by our research group.

An isotropic hardening (IH) model, in which the yield surface is assumed to expand with increasing plastic strain while maintaining a similar shape, is widely used in forming simulations. However, it is well known that the yield surface shape for an actual material changes as work hardening progresses, a phenomenon called differential hardening (DH).^19,20) It is expected that the use of a DH model will improve the accuracy of a forming simulation.²¹⁾ One of the present authors conducted biaxial stress tests on an aluminum alloy sheet^22,23,24) and a steel sheet²⁵⁾ for a large strain range (until the sheet sample fractured) and demonstrated that the analysis accuracy can be improved by applying a DH model with the Yld2000-2d yield function²⁶⁾ to the FEA.

In this study, in addition to IH models based on various yield functions, two types of DH model determined using different parameter identification procedures were created for a hot-rolled steel sheet and the effect of the material model on the prediction accuracy of forming simulations was investigated. The IH and DH models were identified using uniaxial and biaxial tensile test data for the test material. Moreover, a hole expansion forming experiment was conducted to simulate stretch flange forming. The thickness distribution for the test piece after forming was measured. A FEA of hole expansion forming was performed using the created material models. The sheet thickness distribution was calculated and compared with the experimental data. The results show that the DH model effectively improves the accuracy of fracture prediction in hole expansion forming.

2. Experimental Methods

2.1. Test Material

The test material was a hot-rolled steel sheet with a nominal thickness of 1.6 mm. Table 1 shows its mechanical properties. The rolling (RD), transverse (TD), and thickness directions of the material correspond to the x, y, and z axes, respectively. The direction 45° from the RD is referred to as the diagonal direction (DD).

Table 1. Mechanical properties of the test material.

Tensile direction/°	E/GPa	c*/MPa	n*	α*	${\epsilon }_{\text{TS}}^{\text{p}}$**	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~${\epsilon }_{\text{TS}}^{\text{p}}$

**  Logarithmic plastic strain giving the maximum tensile load

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

2.2. Uniaxial Tensile Test

Uniaxial tensile tests were performed at increments of 15° from the RD using JIS13B test pieces. The r-value and uniaxial plastic flow stress measured at θ degrees from the RD are denoted as r_θ and σ_θ, respectively.

2.3. Biaxial Tensile Test

To determine the yield function that accurately reproduces the elastoplastic deformation behavior of the test material in the biaxial stress state, a biaxial tensile test¹¹⁾ with a cruciform test piece and a multiaxial tube expansion test (MTET)¹²⁾ with a circular tube test piece were performed. Figure 1 shows the dimensions of each test piece.

Fig. 1.

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

The geometry of the cruciform test piece and the biaxial tensile test method conformed to ISO16842.¹⁰⁾ A FEA^28,29) confirmed that the stress measurement error can be suppressed to less than 2% by measuring the strain components with strain gauges at the positions shown in the figure.

For the MTET, as-received sheet samples were uniformly bent to form a cylinder and the sheet edges were laser-welded. The RD for type I specimens was in the circumferential direction and that for type II specimens was in the axial direction. See Ref. 12) for details of the test procedures. Because the stress-strain curves measured using the MTET include the effect of prestrain due to bending during the production of the circular tube test piece, they were corrected by the true stress-plastic strain curves measured using a cruciform test piece. For details of the correction method, refer to Ref. 12).

Both cruciform and tubular specimens were subjected to seven linear stress paths, namely σ_x:σ_y = 4:1, 2:1, 4:3, 1:1, 3:4, 1:2, and 1:4, where σ_x and σ_y are the true stress in the RD and TD, respectively. The equivalent plastic strain rate was controlled to be approximately constant at 10⁻⁴ s⁻¹. The MTET was conducted so that the circumferential direction of the tubular test piece was always the maximum principal stress direction. That is, a type I test piece was used for σ_x:σ_y = 4:1, 2:1, and 4:3 and a type II test piece was used for σ_x:σ_y = 3:4, 1:2, and 1:4.

3. Biaxial Stress Test Results

The contours of plastic work in the stress space^19,20) were measured to identify appropriate material models for the test samples subjected to uniaxial and biaxial tension. With the true stress-logarithmic plastic strain curve measured for the RD used as reference data for work hardening, the plastic work per unit volume, W₀, and uniaxial true stress, σ₀, associated with particular values of ${\epsilon }_0^{\text{p}}$ (referred to as the reference plastic strain hereafter) were determined. Next, from the biaxial and uniaxial stress-strain curves, the stress points that give the same plastic work as W₀ were plotted in the principal stress space to determine the contour of plastic work associated with ${\epsilon }_0^{\text{p}}$. Figure 2 shows the measured work contours. We were able to measure the work contours for a maximum value of ${\epsilon }_0^{\text{p}}$ = 0.22.

Fig. 2.

Measured stress points that form the contours of plastic work for various values of ${\epsilon }_0^{\text{p}}$.

To quantitatively evaluate the shape change (DH) of the work contours with ${\epsilon }_0^{\text{p}}$, the nondimensional work contours were determined by dividing the value of the stress points that formed a work contour by the σ₀ belonging to the work contour. Next, we determied the shape ratio, l/l_0.02, for each linear stress path, where l_0.02 is the distance from the origin of the stress space to the stress point (σ_x/σ₀, σ_y/σ₀) that forms the work contour with ${\epsilon }_0^{\text{p}}$ = 0.02, and l is the distance from the origin to the stress point, (σ_x/σ₀, σ_y/σ₀), corresponding to the specified ${\epsilon }_0^{\text{p}}$. If two work contours are similar, l/l_0.02 should be 1 for all stress paths. Figure 3 shows l/l_0.02 versus ${\epsilon }_0^{\text{p}}$. The data for ${\epsilon }_0^{\text{p}}$ ≤ 0.02 are not shown because the shape change of the work contours was unstable due to the influence of the yield elongation of the material. For 0.02 ≤ ${\epsilon }_0^{\text{p}}$ ≤ 0.05, l/l_0.02 changes by up to 3%, and for ${\epsilon }_0^{\text{p}}$ > 0.05, it is almost constant. Thus, the DH for the test sample can be quantitatively evaluated. Such a shape change of work contours cannot be expressed using the IH model.

Fig. 3.

Variation of the shape ratio, l/l_0.02, with ${\epsilon }_0^{\text{p}}$.

4. Material Modeling

IH models were determined by applying the experimental data obtained in the previous section to the von Mises,³⁰⁾ Hill’s quadratic,³¹⁾ and Yld2000-2d²⁶⁾ yield functions. The DH models were determined using the Yld2000-2d yield function. This section describes the derivation for each material model.

4.1. IH Model

4.1.1. Von Mises and Hill’s Yield Functions

The von Mises yield function was identified using σ₀. Hill’s yield function was identified using σ₀ and r-values (r₀, r₄₅, and r₉₀; see Table 1) measured at a nominal strain of 0.10.

4.1.2. Yld2000-2d Yield Function

The Yld2000-2d yield function can be determined using an exponent M and eight parameters α_i (i=1–8). For the IH model based on the Yld2000-2d yield function, a genetic algorithm was used to find the values of M and α_i (i =1–8) that minimize the function f:   

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

The first term is the sum of the difference between the yield surface calculated using the Yld2000-2d yield function and the experimental work contour along each stress path, and the second term is the sum of the difference between the direction of the plastic strain rate, D^p, calculated using the Yld2000-2d yield function and the experimental data for each stress path. l_M,i and l_C,i are the distances along the i^th stress path between the origin of the stress space and the nondimensional work contour and the calculated yield surface, respectively, and β_M,i and β_C,i (units: degrees) are the directions of D^p, measured and calculated using the Yld2000-2d yield function, respectively (see Fig. 4). w_σ,i and w_β,i are weighting factors for each stress path. The values of w_σ,i and w_β,i are shown in Table 2. N (= 14) is the number of stress paths. Seven uniaxial tensile directions (in increments of 15° from the RD) and seven linear stress paths were used for the biaxial tensile test. Regarding the experimental stress points that form a work contour, the flow stress corresponding to ${\epsilon }_0^{\text{p}}$ = 0.22 was used for both uniaxial and biaxial stress points. Regarding the experimental direction of D^p, the r-values at a nominal strain of 0.1 for uniaxial tension and those corresponding to ${\epsilon }_0^{\text{p}}$ = 0.22 for biaxial tension were used.

Fig. 4.

Schematic illustrations for identifying the parameters 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. Comparison of Experimental Data with Values Calculated Using IH Model

Figure 5 compares the measured work contour, directions of D^p, σ_θ, and r_θ with the values calculated using the selected yield functions. The predictions obtained with the Yld2000-2d yield function had the best agreement with the measurements.

Fig. 5.

Measured data for (a) the stress points that form the 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 selected yield functions (IH model). M = 7.3 for the Yld2000-2d.

4.2. DH Model

The DH models were made by approximating the evolution of the M and α_i (i = 1–8) of the Yld2000-2d yield function with plastic strain as a function of ${\epsilon }_0^{\text{p}}$. The following two methods (methods A and B) were used as function approximation methods for the evolution of M and α_i (i = 1–8). Method A, which was developed in Ref. 22), identifies M and α_i (i = 1–8) at the specified ${\epsilon }_0^{\text{p}}$, and approximates each parameter individually as a function of ${\epsilon }_0^{\text{p}}$. In this method, M and α_i (i = 1–8) are optimized discretely for each ${\epsilon }_0^{\text{p}}$, and thus it is difficult to optimize the parameters while reproducing the continuous evolution of M and α_i (i = 1–8) with increasing ${\epsilon }_0^{\text{p}}$. In addition, it takes time to create a model because it is necessary to identify the parameters for each ${\epsilon }_0^{\text{p}}$ each time the weighting coefficients change. Therefore, in this study, we developed method B, in which continuous parameter transitions can be optimized by simultaneously optimizing the approximate expressions for all the specified ${\epsilon }_0^{\text{p}}$.

4.2.1. Creating a DH Model Based on Method A (DH-A)

The procedure for method A was as follows.

i) M and α_i (i = 1–8) were determined at several specified ${\epsilon }_0^{\text{p}}$. The determination procedure is as follows. M was changed in increments of 0.01. For every M, α_i (i = 1–8) were determined using the Newton-Raphson method with eight experimental values (σ₀, σ₄₅, σ₉₀, flow stress σ_b, plastic strain rate ratio $r_{\text{b}}\equiv {\dot{\epsilon }}_y^{\text{p}}/{\dot{\epsilon }}_x^{\text{p}}$ at equibiaxial tension at the specified ${\epsilon }_0^{\text{p}}$, and r₀, r₄₅, and r₉₀ shown in Table 1). The error between the identified yield locus and the experimental work contour was evaluated using Eq. (1) to determine M that minimizes the function f. The weighting factors were w_σ,i = 1 and w_β,i = 0.3.

ii) The change in M obtained in i) was approximated as a function of ${\epsilon }_0^{\text{p}}$:   

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

iii) Using the M approximated by Eq. (2), the α_i values (i = 1–8) for each ${\epsilon }_0^{\text{p}}$ were determined using a genetic algorithm (the calculation procedure was the same as that described in Section 4.1.2).

iv) The changes in α_i (i = 1–8) obtained in iii) were approximated using the following equation as a function of ${\epsilon }_0^{\text{p}}$:   

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

See Table A1 for the values of the coefficients A, B, C, and D in Eqs. (2) and (3).

4.2.2. Creating a DH Model Based on Method B (DH-B)

The procedure for method B was as follows.

i) M and α_i (i = 1–8) were approximated as a function of ${\epsilon }_0^{\text{p}}$ using the following equations (the same procedure as that in Section 4.2.1).   

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

ii) With the parameters A, B, C, and D in Eq. (4) as initial values, the objective function f_B of Eq. (5) was minimized using a genetic algorithm.   

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

See Table A2 for the values of parameters A, B, C, and D.

4.2.3. Comparison of Experimental Data with Values Calculated Using DH Model

Comparisons of the work contours at ${\epsilon }_0^{\text{p}}$ = 0.005 and 0.22, the directions of D^p, σ_θ, and r_θ with those calculated using the DH model are shown in Fig. 6. For reference, the calculated values using the IH model based on the Yld2000-2d yield function are also shown in the figure.

Fig. 6.

Measured data for (a) the stress points that form the 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).

Because the IH model was based on the experimental data at ${\epsilon }_0^{\text{p}}$ = 0.22, there is a discrepancy with the experimental data at ${\epsilon }_0^{\text{p}}$ = 0.005. The DH model accurately reproduces the evolution of the work contours and σ_θ with increasing strain. Regarding r_θ, although there is a slight deviation from the experimental data for both DH-A and DH-B, the tendency of the r-value reaching its maximum at the DD is qualitatively reproduced.

5. Hole Expansion Experiment and Finite Element Analysis

The usefulness of the material models created in the previous section was investigated by applying them to a hole expansion simulation and comparing the predicted thickness strain distribution with the measured one.

5.1. Hole Expansion Experiment

Figure 7(a) shows a cross-sectional view of the experimental apparatus used for hole expansion forming. A circular blank with a diameter of 215 mm was used. A hole with a diameter of 30 mm was made in the center of the blank using wire electric discharge machining. The punch diameter was 100 mm and the radius of the punch and die profile was 15 mm. The punch stroke speed was 0.1 mm/s. A grid with 10° increments in the circumferential direction and 2-mm increments in the radial direction was scribed on the blank and used to identify the initial coordinates for thickness measurements after the test. A Teflon sheet coated with Vaseline on both sides was inserted between the punch and the blank for lubrication. The punch stroke, h, was 28 mm.

Fig. 7.

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

5.2. FEA of Hole Expansion Forming

For the FEA of the hole expansion forming, the static implicit software ABAQUS/Standard 2017 was used. An outline of the FEA model is shown in Fig. 7(b). The die was an analytical rigid body. Four-node reduced integration shell elements (S4R) were used for the blank. The number of integration points in the thickness direction was seven. One quarter of the blank was analyzed considering symmetry. The blank was divided into finite elements every 5° in the circumferential direction and every 1 mm in the radial direction. The initial thickness of the blank was 1.574 mm, which was the average initial sheet thickness of the test sample used in the experiment, and the initial hole diameter was 30 mm. Under the assumption of no material inflow from the apex of the bead, the node displacement on the outer edge of the blank with a radius of 97.5 mm was fixed. The friction coefficient was 0.03 at the interface between the blank and the punch and 0.2 at the interface between the die and the blank. Swift’s hardening equation was used as the work hardening equation based on the RD uniaxial tensile test results (see Table 1).

5.3. Comparison of Calculated and Experimental Results

5.3.1. IH Model

Figure 8 shows a comparison of the calculated and experimental values of the plastic thickness strain, ${\epsilon }_z^{\text{p}}$, in the circumferential direction at a position of 2 mm distant from the hole edge on the initial blank. The experimental values of ${\epsilon }_z^{\text{p}}$ exhibited a maximum for the DD, a local minimum for the RD and TD, and an absolute minimum for the RD. The von Mises yield function, which assumes isotropy, naturally shows a uniform thickness reduction. For Hill’s yield function, the phases of the maximum and minimum values are opposite to those for the experimental values; the accuracy of the material model is thus low. In contrast, for the Yld2000-2d yield function, the positions of the extrema are in agreement with the experimental values. However, the minimum value is predicted to be at the TD; therefore, the experimental result of ${\epsilon }_z^{\text{p}}$ reaching its absolute minimum at the RD is not reproduced.

Fig. 8.

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

Figure 9 compares the experimental values of the RD, DD, and TD distributions of ${\epsilon }_z^{\text{p}}$ with those calculated using the IH models. The Yld2000-2d yield function has a maximum difference of $\left|\mathrm{\Delta }{\epsilon }_z^{\text{p}}\right|$ ≈ 0.05 from the experimental value in the DD; the difference from those in the RD and TD is within $\left|\mathrm{\Delta }{\epsilon }_z^{\text{p}}\right|$ = 0.02, accurately reproducing the experimental data. In addition, the tendency of the minimum value being at a position 4 to 6 mm away from the hole edge in the RD and DD is well reproduced.

Fig. 9.

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

5.3.2. DH Model

Figure 10 compares the experimental values of the circumferential ${\epsilon }_z^{\text{p}}$ distribution with those calculated using the DH models. Both DH-A and -B predict the absolute minimum in the RD and the local minimum in the TD; the predicted values are almost the same as the experimental values. Therefore, the analysis accuracy is improved compared to that obtained with the IH model. Both DH-A and -B models underestimate the amplitude of ${\epsilon }_z^{\text{p}}$. In the DD, the reduction of the sheet thickness is overestimated by about Δ${\epsilon }_z^{\text{p}}$ ≈ 0.05 to 0.04. The result for DH-B is slightly closer to the experimental value.

Fig. 10.

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

Figure 11 compares the experimental values of the radial ${\epsilon }_z^{\text{p}}$ distribution with those calculated using the DH models. In the RD, DD, and TD, the difference between the three material models (IH, DH-A, and DH-B) is Δ${\epsilon }_z^{\text{p}}$ ≈ 0.01, which is approximately on the same order as the variation in the experiment. In the RD and TD, both DH models reproduce the experimental data with an error in the range of Δ${\epsilon }_z^{\text{p}}$ ≈ 0.005 to 0.03. The results for DH-A are the closest to the experimental values in the RD, and those for IH and DH-B are the closest in the TD. In the DD, the results for IH are the closest to the experimental values, but all material models overestimate the sheet thickness reduction in the range 0.02 ≤ Δ${\epsilon }_z^{\text{p}}$ ≤ 0.06, and the deviation from the experimental value is larger than those in the RD and TD. These observations do not conclusively indicate which material model is best in terms of the accuracy of predicting the radial ${\epsilon }_z^{\text{p}}$ distribution.

Fig. 11.

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

5.4. Discussion

The effect of material models on the results and accuracy of the hole expansion forming simulation are summarized as follows. In the comparison between IH models, the Yld2000-2d yield function most accurately reproduced the experimental values of the ${\epsilon }_z^{\text{p}}$ distribution. Furthermore, using the model that reproduces the DH of the test sample, we were able to accurately identify the minimum position of ${\epsilon }_z^{\text{p}}$ at the hole edge. This is thought to be due to the DH model accurately reproducing the evolution of the circumferential distribution of σ_θ with ${\epsilon }_0^{\text{p}}$. Figure 12 compares the experimental values of σ_θ with those calculated using DH-B. In the low strain range of 0.005 ≤ ${\epsilon }_0^{\text{p}}$ ≤ 0.01, σ_θ monotonically increases from the RD to the TD; however, in the range of 0.02 ≤ ${\epsilon }_0^{\text{p}}$, σ_θ becomes minimum in the DD; DH-B can accurately reproduce these experimental data. Because the hole edge is in a uniaxial tensile stress state, it is presumed that the DH model, which has excellent reproduction accuracy of the development of σ_θ, improved the prediction accuracy of the development of ${\epsilon }_z^{\text{p}}$ near the hole edge compared to that of the IH model. When the punch was further raised in the hole expansion test, the blank fractured at a position inside the hole edge in the RD at h = 32 mm, as shown in Fig. 13. From this, it is considered that the DH model that accurately reproduces the minimum value of ${\epsilon }_z^{\text{p}}$ in the RD is effective in predicting the fracture of stretch flange forming. This result is consistent with the findings reported in Refs. 18) and 24).

Fig. 12.

Variation of σ_θ/σ₀ with ${\epsilon }_0^{\text{p}}$ compared with that calculated using the DH-B model.

Fig. 13.

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

Even with the DH model, the thickness reduction in the DD was overestimated in the range of 0.02 ≤ Δ${\epsilon }_z^{\text{p}}$ ≤ 0.06, and the deviation from the experimental values was larger than that in the RD and TD (Fig. 11(b)). This means that the material model was not perfect. In this study, the parameters for the Yld2000-2d yield function were determined to accurately reproduce the work contour shape in the σ_x−σ_y plane stress space, the directions of D^p, and the distributions of r_θ and σ_θ, which represent the deformation and stress states near the hole edge. To establish a more accurate material model, it is necessary for the material model to reproduce the three-dimensional yield surface shape and the development of D^p in the σ_x – σ_y − σ_xy space with high accuracy. To this end, additional material tests, described below, should be carried out. Because the anisotropy of the test material is relatively small, let us assume that the test sample is isotropic to simplify the discussion. Then, the stress state along the DD during hole expansion forming will be on the σ_x = σ_y plane in the σ_x – σ_y − σ_xy stress space. To reproduce these stress states, biaxial tensile tests should be performed with the arm direction of the cruciform test piece parallel to the DD. If a material model that can reproduce the plastic deformation characteristics along the DD obtained in this way can be made, the calculated value of the thickness distribution in the DD will be closer to the experimental data. Another method for measuring the evolution of the yield surface in the σ_x – σ_y − σ_xy stress space is to perform plane strain tensile tests with the principal stress axis in several directions to represent the deformation modes at positions some distance from the hole edge, where the principal strain rate ratio is ${\dot{\epsilon }}_1^{\text{p}}\text{:}{\dot{\epsilon }}_2^{\text{p}}$ = 1:0.³²⁾

6. Conclusions

Biaxial stress tests were performed on a 440-MPa hot-rolled steel sheet and the DH behavior was precisely measured. In addition to the IH model, two types of DH model that reproduce the DH behavior were created. A hole expansion forming simulation was performed using the three types of material model and the effect of the material model on the accuracy of the forming simulation was evaluated by comparing the calculated values with the experimental values. The following findings were obtained.

(1) The DH behavior of the test material was quantitatively evaluated from the work contours in the first quadrant of the σ_x − σ_y stress space (σ_z = 0) and the directions of D^p measured using biaxial stress tests with seven linear stress paths, and from uniaxial tensile test data taken at 15° increments in the tensile direction.

(2) DH models for the test sample were created by expressing the evolution of the exponent M and α_i (i = 1–8) of the Yld2000-2d yield function as a function of ${\epsilon }_0^{\text{p}}$.

(3) With the developed DH models in the hole expansion forming simulation, we were able to accurately reproduce the experimental result of ${\epsilon }_z^{\text{p}}$ becoming a minimum in the RD. Therefore, the DH model is effective for improving the accuracy of fracture prediction in stretch flange forming.

(4) To further improve the accuracy of the material model, it is necessary to reproduce the three-dimensional yield surface shape and the development of D^p in the σ_x − σ_y − σ_xy stress space with high accuracy. For that purpose, the following material tests should be carried out and a material model that can reproduce the obtained plastic deformation characteristics should be constructed: (i) biaxial tensile tests with the arm direction of the cruciform test piece parallel to the DD and various stress ratios; (ii) plane strain tensile tests (principal strain rate ratio of ${\dot{\epsilon }}_1^{\text{p}}\text{:}{\dot{\epsilon }}_2^{\text{p}}$ = 1:0) with various principal stress directions.

Acknowledgement

We would like to express our deepest gratitude to Dr. Nobuyasu Noman (Unipres Corporation) for providing the test material.

References

• 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), 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,  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 et al., 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.  Ogihara 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: Nachr. Ges. Wiss. Göttingen, Math.-Phys. Kl., (1913), 582. http://eudml.org/doc/58894, (accessed 2020-03-13).
• 31)   R.  Hill: Proc. R. 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

Appendix

Table A1 shows the parameters in Eqs. (2) and (3) used for the DH-A model developed in Section 4, and Table A2 shows the parameters in Eq. (4) used for the DH-B model. The variation in α₅ and α₇ for DH-A with the evolution of strain was small, and thus the average value was taken.

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\left({\epsilon }_0^{\text{p}}\right)=\left(A-B\right)/\left[1+\text{exp}\left\{\left({\epsilon }_0^{\text{p}}-C\right)/D\right\}\right]+B$ (Eq. (2))

**  Approximated using ${\alpha }_i\left({\epsilon }_0^{\text{p}}\right)=A-B\text{exp}\left(-C{\epsilon }_0^{\text{p}}\right)-\left\{D/\left(0.001+{\epsilon }_0^{\text{p}}\right)\right\}$ (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,\mathrm{\hspace{0.17em} } {\alpha }_i(i=1-8)=A{\epsilon }_0^{\text{p}}\left({\epsilon }_0^{\text{p}}-B\right)\text{exp}\left(C{\epsilon }_0^{\text{p}}\right)+D$ (Eq. (4))