Introduction

Background

The sheet metal forming (SMF) industry is highly competitive and imposes heavy demands regarding cost, lead time, and quality, et.c. SMF is utilized in fields such as vehicles, energy, aerospace, and home appliances. In these industries, finite element (FE) simulations are utilized as a valuable tool during the development of new products that are manufactured with stamping dies. Such simulations have been key to the rapid development in the industry in the last decades. Thorough background information on sheet metal forming processes can be found in [1], which focuses on theoretical modelling of the forming process and material models.

An important part of reliable sheet metal forming simulations is the theoretical description of the thin metal sheet, where hardening and yield criteria are significant [1]. The focus of this paper is the mathematical description, and calibration, of the yield criteria. The most commonly applied and widely known yield criteria are Tresca and von Mises for isotropic material behavior. However, many metals exhibit anisotropic behavior due to manufacturing processes and require a mathematical approximation to reliably describe their behavior. The first yield criterion for anisotropic materials was proposed by von Mises [1,2,3,4,5,6] and later modified by Hill in 1948 [6] to describe orthotropic materials. Hill (1948) as well as Tresca and von Mises are still widely referenced within both academia and the industry. The benefits of these yield criteria are that they are relatively easy to understand and require only a few experiments to determine the yield surface parameters. However, a disadvantage is that they cannot describe material anisotropy for new and advanced materials with sufficient accuracy for the high demands of simulations in today’s industry and academia.

Researchers have developed several modifications of the previously mentioned yield criteria as well as completely new models. Currently, many advanced yield criteria exist that were specifically developed for plane stress assumptions and are often used in sheet metal forming simulations (e.g. Barlat 2000 [7, 8], Vegter [9, 10], and Banabic-Balan-Comsa [1, 11, 12]). Several modifications have been made to most of these advanced yield criteria for thin metal sheets, and even though they derive from different theories and backgrounds, they are equivalent to each other under certain conditions [13].

Experimental input are required to determine the coefficients of all the yield criteria. Typical input parameters are yield stresses for different stress conditions and coefficients of anisotropy. The yield criteria that were applied in this research is BBC05 [1], which are commonly used in the sheet metal forming industry and implemented in several FE software. Up to eight parameters can be employed to determine this yield surface, including three yield stresses in three directions of the metal sheet as well as the biaxial stress. The other four coefficients are coefficients of anisotropy that are often called r values [1].

The yield surface exponent: k

An exponent, k, is also present in the expression of the BBC yield surface. This exponent is not determined directly from experiments but traditionally set to 3 for steel alloys (BCC) and 4 for aluminum alloys (FCC) in accordance with the crystallographic structure of the material. These values for the exponent should only be taken as guiding principles. It is visualized in [14] that the exponent can vary a lot for metallic materials. The utilized material model in [14] is Barlat 2000 [7, 8], equivalent to BBC05 [15]. Practical experience at Volvo Cars also supports exponent values other than 3 and 4. The variation of the exponent k is also discussed in [15] along with a comparison between the Vegter yield critera and models such as Barlat 2000, and Hill (1948). The assumption that k should be 3 or 4 was originally made for an extension of the isotropic von Mises quadratic yield function [16]. Since both the mathematical formulations differ per yield criteria, and most metal sheet materials are not isotropic, the values of 3 and 4 for k should not be accepted as absolute truth. For an accurate description of many materials, this coefficient is required to adopt values significantly higher or lower than the proposed values of 3 and 4. Figure 1 depicts a yield surface with the same experimental input for stresses and r values but different values for the exponent k.

Fig. 1
figure 1

First quadrant of yield surfaces with different exponent k; values of exponent correspond to results in later sections

Full size image

An example of the influence of the exponent k in simulations of industrial parts is depicted in Fig. 2. The part is an XC90 front inner door from the running production at the Volvo Cars plant in Olofström. The simulations was carried out in accordance with the methods presented in [17], and the only difference between the depicted results are the value of k. The figure demonstrate the amount of over- and under prediction of the major strains in the simulated door. The simulated strains are compared to experimental values according to the following equation:

$$ Strain\ prediction= Simulated\ strains- Measured\ strains $$
(1)
Fig. 2
figure 2

Simulation of XC90 front inner door with k=3.5 and 2.6, deviation in strain according to Eq. 1

Full size image

The simulations reveal similar strain predictions; however, in sensitive areas such as vertical walls, the variation of k has a relatively large influence. It is imperative to ensure reliability of the strain predictions in these, and other sensitive areas, since cracks often form here in running production and during tryout.

Research focus and purpose

The first purpose of this paper is to present the BBC05 yield criterion with non-integer exponent. The BBC05 yield criterion is accompanied by a description of how to implement it in an explicit FE code. This description of the numerical implementation can also be used by other engineers and researchers as a basis for their own implementation of BBC05 or other yield criteria. There are many papers presenting material research using user sub routines, however, descriptions of the numerical implementation in FE codes is more seldom described.

The second purpose is to establish a robust and accurate procedure for determining the value of the yield surface exponent k. This procedure is carried out through inverse modeling of nakajima experiments, also used for FLC evaluation. The utilized nakajima samples corresponds to equi-biaxial strain and plane strain. Plane strain experiments have previously been used in the literature for yield surface calibration, e.g. [18]. Yield surface calibration with Nakajima tests has also been studied in [19, 20] where topics such as flow rules, crystal plasticity, yield loci definitions among other topics are also included.

All sheet metal forming simulations at Volvo Cars apply the BBC05 yield criteria with non-integer exponent instead of the BBC05 version with integer exponent in [1]. Isotropic hardening is also used at Volvo Cars and throughout this paper. The use of material models with non-integer exponents introduces more flexibility to the material model during the calibration of the yield surface. This paper presents a modification of the BBC05 with non-integer exponent.

Volvo Cars conducts tensile tests in three directions to determine hardening curves and yield stresses, together with viscous bulge tests [21] to also determine biaxial data for the BBC05 yield criteria. For determining k and fine tune the material model, additional experimental data need to be evaluated as described in the section “Calibration of yield surface” in this paper. In Nakajima tests, the strain and punch force is sensitive to the value of k. Therefore, it can provide data for determining the exponent k in the BBC05 model. However, friction between sheet and tools also influences the Nakajima tests, which can cause ambiguities in the determination of the exponent k. The influence of friction level and the value of k similarly influence the strains and forces, this will be exemplified by simulations of Nakajima tests without lubrication (often called LDH test [1]). Consequently, it can be difficult to separate the variables friction and k unless the experiments are performed to ensure that the friction level is known, or so low that it can be neglected.

BBC05 yield surface with integer exponent

The following section summarizes the most important equations and concepts for explaining the BBC05 material model. A full mathematical background together with identification procedures for the BBC05 yield surface can be found in [1]. The following conventions are used: Greek indices take the values 1 and 2, while Latin indices take the values 1, 2, and 3. The directions 1, 2, and 3 in the selected basis system for strains and stresses respectively correspond to rolling direction, transverse direction, and normal direction of the thin sheet metal.

Flow rule

The following equation describes the yield surface for a sheet metal that behaves as a plastically orthotropic membrane:

$$ \Phi \left({\sigma}_{\alpha \beta},Y\right)=\overline{\sigma}\left({\sigma}_{\alpha \beta}\right)-Y=0 $$
(2)

where \( \overline{\sigma}>0 \) is the BBC05 equivalent stress, Y > 0 is a yield parameter that can be chosen without any constraint, and σαβ = σβα are planar components of the stress tensor. Restrictions that arise from plane stress conditions are

$$ {\sigma}_{3i}={\sigma}_{i3}=0. $$
(3)

The flow rule is expressed as

$$ {\overset{.}{\varepsilon}}_{\alpha \beta}^p=\overset{.}{\lambda}\frac{\mathrm{\partial \Phi }}{\partial {\sigma}_{\alpha \beta}} $$
(4)

which is an associated flow rule, and \( \overset{.}{\lambda}\ge 0 \) is a plastic multiplier. According to the consistency condition, the stress point always remains on the yield surface during plastic loading. For yield surfaces with an associated flow rule that obey the consistency condition the following relationship holds

$$ \overset{.}{\lambda }=\overset{.}{\overline{\varepsilon^p}} $$
(5)

where \( \overset{.}{\overline{\varepsilon^p}} \) is the equivalent plastic strain rate. This relationship is important when numerically implementing the BBC05 model.

Restrictions for the out-of-plane components of the plastic strain-rate tensor are as follows:

$$ {\overset{.}{\varepsilon}}_{3\alpha}^p={\overset{.}{\varepsilon}}_{\alpha 3}^p=0 $$
(6)
$$ {\overset{.}{\varepsilon}}_{33}^p=-\left({\overset{.}{\varepsilon}}_{11}^p+{\overset{.}{\varepsilon}}_{22}^p\right) $$
(7)

Equivalent stress

The equivalent stress for BBC05, \( \overline{\sigma} \), is expressed as

$$ \overline{\sigma}={\left[a{\left(\Lambda +\Gamma \right)}^{2k}+a{\left(\Lambda -\Gamma \right)}^{2k}+b{\left(\Lambda +\Psi \right)}^{2k}+b{\left(\Lambda -\Psi \right)}^{2k}\right]}^{\frac{1}{2k}} $$
(8)

where a, b > 0, and k is a positive integer. In addition, Γ, Λ, and Ψ are functions that depend on the planar components of the stress tensor, and a, b, L, M, N, P, Q, R, and k are all material parameters.

$$ {\displaystyle \begin{array}{c}\kern1.4em \varGamma =L{\sigma}_{11}+M{\sigma}_{22}\\ {}\varLambda =\sqrt{{\left(N{\sigma}_{11}-P{\sigma}_{22}\right)}^2+{\sigma}_{12}{\sigma}_{21}}\\ {}\varPsi =\sqrt{{\left(Q{\sigma}_{11}-R{\sigma}_{22}\right)}^2+{\sigma}_{12}{\sigma}_{21}}\end{array}} $$
(9)

Application of the flow rule that was given for BBC05 requires the partial derivatives of the function Φ with respect to the planar components of the stress tensor, which is expressed by the following equation:

$$ {\displaystyle \begin{array}{c}\frac{\partial \varPhi }{\partial {\sigma}_{\alpha \beta}}=\frac{\partial \overline{\sigma}}{\partial {\sigma}_{\alpha \beta}}\\ {}\frac{\partial \varPhi }{\partial {\sigma}_{\alpha \beta}}=\frac{\partial \overline{\sigma}}{\partial \varGamma}\frac{\partial \varGamma }{\partial {\sigma}_{\alpha \beta}}+\frac{\partial \overline{\sigma}}{\partial \varLambda}\frac{\partial \varLambda }{\partial {\sigma}_{\alpha \beta}}+\frac{\partial \overline{\sigma}}{\partial \varPsi}\frac{\partial \varPsi }{\partial {\sigma}_{\alpha \beta}}\end{array}} $$
(10)

Equation 8 together with Eq. 10 gives

$$ \frac{\partial \overline{\sigma}}{\mathrm{\partial \Gamma }}=\frac{a}{\overline{\sigma^{2k-1}}}\left[{\left(\Lambda +\Gamma \right)}^{2k-1}-{\left(\Lambda -\Gamma \right)}^{2k-1}\right] $$
$$ \frac{\partial \overline{\sigma}}{\mathrm{\partial \Lambda }}=\frac{1}{\overline{\sigma^{2k-1}}}\left\{a\left[{\left(\Lambda +\Gamma \right)}^{2k-1}+{\left(\Lambda -\Gamma \right)}^{2k-1}\right]+b\left[{\left(\Lambda +\Psi \right)}^{2k-1}+{\left(\Lambda -\Psi \right)}^{2k-1}\right]\right\} $$
$$ \frac{\partial \overline{\sigma}}{\mathrm{\partial \Psi }}=\frac{b}{\overline{\sigma^{2k-1}}}\left[{\left(\Lambda +\Psi \right)}^{2k-1}+{\left(\Lambda -\Psi \right)}^{2k-1}\right] $$
(11)

BBC05 material model with non-integer exponent

As discussed in sub-section The yield surface exponent: k, the material parameter k is an integer that has been traditionally selected in accordance with the crystallographic structure of the material. k = 3 for BCC material, while k = 4 for FCC materials. Therefore, k is held as a constant throughout the identification process for the constants of the BBC05 material model. However, practical experience with inverse modelling of a large number of materials at Volvo Cars revealed that the BBC05 model can often describe the behavior of a material with significantly higher precision if k is allowed to be a real positive number instead of a positive integer.

The BBC05 material model can be expressed with non-integer exponent by calculating the effective stress with moduli instead of plain parentheses according to Eq. 12.

$$ \overline{\sigma}={\left[a{\left|\Lambda +\Gamma \right|}^{2k}+a{\left|\Lambda -\Gamma \right|}^{2k}+b{\left|\Lambda +\Psi \right|}^{2k}+b{\left|\Lambda -\Psi \right|}^{2k}\right]}^{\frac{1}{2k}} $$
(12)

The identification procedure of the constants can be performed in the same way as for BBC05 with integer exponent as described in [1], modified equations for BBC05 with non-integer exponent can be found in section 12 Appendix of this paper. Differences are inevitable in the expressions for the partial derivatives of \( \overline{\sigma} \) with respect to Λ, Γ, and Ψ. In the BBC05 version with non-integer exponent, Eq. 11 becomes:

$$ {\displaystyle \begin{array}{c}\frac{\partial \overline{\sigma}}{\partial \varGamma }=\frac{a}{\overline{\sigma^{2k-1}}}\left[{\left\langle \varLambda +\varGamma \right\rangle}^{2k-1}-{\left\langle \varLambda -\varGamma \right\rangle}^{2k-1}\right]\\ {}\frac{\partial \overline{\sigma}}{\partial \varLambda }=\frac{1}{\overline{\sigma^{2k-1}}}\left\{a\left[{\left\langle \varLambda +\varGamma \right\rangle}^{2k-1}+{\left\langle \varLambda +\varGamma \right\rangle}^{2k-1}\right]+b\left[{\left\langle \varLambda +\varPsi \right\rangle}^{2k-1}+{\left\langle \varLambda -\varPsi \right\rangle}^{2k-1}\right]\right\}\\ {}\frac{\partial \overline{\sigma}}{\partial \varPsi }=\frac{b}{\sigma^{2k-1}}\left[{\left\langle \varLambda +\varPsi \right\rangle}^{2k-1}+{\left\langle \varLambda -\varPsi \right\rangle}^{2k-1}\right]\end{array}} $$
(13)

where the following notation is used

$$ {\left\langle x\right\rangle}^y={\left|x\right|}^y\mathit{\operatorname{sgn}}(x) $$
(14)

Implementation of BBC05 as a subroutine in LS Dyna

The BBC05 material model can be implemented in any FE software. This section describes its implementation in LS Dyna [22] for a solver that uses explicit time integration. However, the outlined procedure is quite general and can be modified for other software or time integration procedures. It can also be used for any other yield surface with an associated flow rule and plane stress conditions. Here, plane stress is assumed and the integration of the constitutive equation is carried out by the radial return algorithm. The reader is referenced to [1, 23] for a deeper understanding of yield surfaces, non-linear FE method and different methods for integration of constitutive equations. User input for the material subroutine are the plastic hardening curve with plastic and elastic material parameters. Means for determining the material coefficients from experimental values have been specified in [1]. The subroutine in LS Dyna functions as follows:

Strain increments, ∆εαβ, and current plastic strain, \( \overline{\varepsilon^p} \), are passed to the subroutine from the FE solver, together with stresses, σαβ, at end of last increment. The plastic strain, \( \overline{\varepsilon^p} \), is used to find the current yield stress, Y, and the slope of the plastic hardening curve, Ep, through the LS Dyna function crvval. (Note: older LS-Dyna versions pass in the strain rate, \( {\overset{.}{\varepsilon}}_{\alpha \beta} \), instead. Users should always be careful and modify equations accordingly depending on the variable and solver or software that are employed).

Trial stresses at the end of the time increment,\( {\sigma}_{\alpha \beta}^{tr} \), in the plane of the metal sheet are then calculated with Hooke’s law for linear elastic plane stress according to

$$ {\displaystyle \begin{array}{c}{\sigma}_{11}^{tr}={\sigma}_{11}^t+{E}_n\left(\Delta {\varepsilon}_{11}+v\Delta {\varepsilon}_{22}\right)\\ {}{\sigma}_{22}^{tr}={\sigma}_{22}^t+{E}_n\left({\Delta \varepsilon}_{22}+v{\Delta \varepsilon}_{11}\right)\\ {}{\sigma}_{12}^{tr}={\sigma}_{12}^t+G\Delta {\varepsilon}_{12}\end{array}} $$
(15)

and the BBC05 effective stress based on the trial stresses, \( {\overline{\sigma}}^{tr} \), are calculated according to Eq. 10. Superscript t indicates the beginning of the current time increment. v is Poisson’s ratio, G is the Shear modulus, and En is the modified Young’s modulus according to

$$ G=\frac{E}{2\left(1+v\right)} $$
(16)
$$ {E}_n=\frac{E}{1-{v}^2} $$
(17)

Equation 2 then incorporates the in-plane stresses to check if the material is yielding.

In Alternative 1, the material does not yield: \( {\overline{\sigma}}^{tr}\le Y \). In this case, the new stresses are updated according to

$$ {\displaystyle \begin{array}{c}{\sigma}_{11}^{t+\Delta t}={\sigma}_{11}^{tr}\\ {}{\sigma}_{22}^{t+\Delta t}={\sigma}_{22}^{tr}\\ {}{\sigma}_{12}^{t+\Delta t}={\sigma}_{12}^{tr}\end{array}} $$
(18)

Where superscript t + ∆t the end of the current time increment. The out-of-plane strain increment is set to

$$ {\Delta \varepsilon}_{33}=-\frac{v}{E}\left(\Delta {\sigma}_{11}+\Delta {\sigma}_{22}\right) $$
(19)

where \( \Delta {\sigma}_{\alpha \beta}={\sigma}_{\alpha \beta}^{t+\Delta t}-{\sigma}_{\alpha \beta}^t \)

Alternative 2 is that the material is yielding: \( {\overline{\sigma}}^{tr}>Y \). The following system of equations must then hold true for elasto-plasticity. Note: All the partial derivatives, \( \frac{\partial \overline{\sigma}}{\partial {\sigma}_{\alpha \beta}} \), in Eq. 20 are evaluated at the end of the current time increment. This is indicated by writing \( {\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{\alpha \beta}}\right|}_{t+\Delta t} \) and is important if larger time steps are used such as in implicit simulations, it works for both explicit and implicit time integrations.

$$ \left\{\begin{array}{c}\Delta {\sigma}_{11}={E}_n\left(\Delta {\varepsilon}_{11}-\Delta \overline{\varepsilon^p}{\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{11}}\right|}_{t+\Delta t}\right)+v{E}_n\left(\Delta {\varepsilon}_{22}-\Delta \overline{\varepsilon^p}{\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{22}}\right|}_{t+\Delta t}\right)\\ {}\Delta {\sigma}_{22}=v{E}_n\left(\Delta {\varepsilon}_{11}-\Delta \overline{\varepsilon^p}{\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{11}}\right|}_{t+\Delta t}\right)+{E}_n\left(\Delta {\varepsilon}_{22}-\Delta \overline{\varepsilon^p}{\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{22}}\right|}_{t+\Delta t}\right)\\ {}\Delta {\sigma}_{12}=G\left(\Delta {\varepsilon}_{12}-2\Delta \overline{\varepsilon^p}{\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{12}}\right|}_{t+\Delta t}\right)\end{array}\right. $$
(20)

Where the plastic strain increment, \( \overline{\varepsilon^p} \), in Eq. 20 is calculated as

$$ \Delta \overline{\varepsilon^p}=\frac{{\overline{\sigma}}^{tr}-Y}{E^p} $$
(21)

The system of equations in Eq. 20 is solved by finding a solution with the trial stresses as variables. The Newton-Raphson or secant method is commonly applied with the elastic trial stresses as a starting point.

When the solution for the trial stresses converges, the yield stress is updated according to

$$ Y={\overline{\sigma}}^{tr} $$
(22)

Where \( {\overline{\sigma}}^{tr} \) is calculated with the trial stresses from after the solution for the equation system in Eq. 20 have converged. The in-plane stresses are updated by

$$ {\displaystyle \begin{array}{c}{\sigma}_{11}^{t+\Delta t}={\sigma}_{11}^{tr}\\ {}{\sigma}_{22}^{t+\Delta t}={\sigma}_{11}^{tr}\\ {}{\sigma}_{12}^{t+\Delta t}={\sigma}_{12}^{tr}\end{array}} $$
(23)

In addition, the plastic strain is updated with

$$ \overline{\varepsilon^p}=\overline{\varepsilon^p}+\Delta \overline{\varepsilon^p} $$
(24)

The element thickness strain increment is calculated as follows:

$$ \Delta {\varepsilon}_{33}=-\left(\Delta \overline{\varepsilon^p}{\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{11}}\right|}_{t+\Delta t}+\Delta \overline{\varepsilon^p}{\left.\frac{\partial \overline{\sigma}}{\partial {\sigma}_{22}}\right|}_{t+\Delta t}\right)-\frac{v}{E}\left(\Delta {\sigma}_{11}+\Delta {\sigma}_{22}\right) $$
(25)

Remark on out-of-plane stresses in LS-Dyna: Under plane stress assumptions, the out-of-plane stresses should be 0; however, σ13 and σ23 must be updated in LS Dyna, which occurs elastically in this case. The out-of-plane stresses are updated according to

$$ {\displaystyle \begin{array}{c}{\sigma}_{33}=0\\ {}{\sigma}_{13}={\sigma}_{13}+2\upkappa \mathrm{G}\Delta {\upvarepsilon}_{13}\\ {}{\sigma}_{23}={\sigma}_{23}+2\upkappa \mathrm{G}\Delta {\upvarepsilon}_{23}\end{array}} $$
(26)

where κ is a shear correction factor for shell elements depending on the element formulation, the recommended value is often \( \frac{5}{6} \). An overview of the BBC05 implementation in LS-Dyna is visualized in Fig. 3, verification of the implementation is discussed in sub section Model differences and verification of BBC05 UMAT implementation.

Fig. 3
figure 3

Overview of material subroutine in LS-Dyna for BBC05 implementation

Full size image

BBC05 experimental input

For calibration of the BBC05 material model, Volvo Cars relies on both direct experimental input and inverse modeling. The following tests were performed at RISE IVF, Olofström, Sweden.

The material presented in this paper is a hot-dip galvanized steel of grade CR4 with a thickness of 0.71 mm and Fuchs RP Anticorit 4107S lubrication. The measured material parameters are as follows:

Tensile tests

Tensile tests, performed according to standard EN 10002–1 [24] with an ARAMIS digital image correlation (DIC) system, determine the yield stresses and r values in the rolling direction, the transverse direction, and 45° from the rolling direction by evaluation in the ARAMIS software. The tests also directly establish the first part of the plastic hardening curve up to ultimate tensile stress. Tests are evaluated with the publically available ARAMIS script “Tensile Test v630” from GOM.(Table 1)

Table 1 Mechanical material parameters for sheet material used in simulation
Full size table

Bulge test

A viscous bulge test is used to calculate the biaxial yield stress, and first part of the extension of the plastic hardening curve after diffuse necking. Experimental procedure, test equipment, and evaluation methods at Volvo Cars are described in [21]. The first part of a Bulge curve should not be used since the evaluation of the results are based on the curvature of the specimen. The curvature is very low in the beginning of the Bulge test and therefore distorts the results. The beginning of the Bulge curve is therefore not presented in Fig. 4. The Tensile curve and the Bulge curve are joined in the area where they overlap. The biaxial r-value is evaluated from the strains at the pole of the bulge according to Eq. 27. The final plastic hardening curve is then obtained by weighting a Hockett-Sherby extension of the tensile curve together with a Voce extension. The Hockett-Sherby extension receives a weight of 0.9, while the Voce extension is assigned a weight of 0.1, which produces an extension of the tensile test curve that follows the curve that is given by the bulge test.

$$ {r}_{biax}=\frac{\varepsilon_{22}}{\varepsilon_{11}} $$
(27)
Fig. 4
figure 4

Tensile curve and transformed bulge curve

Full size image

Nakajima tests (FLC and LDH)

The Volvo Cars Nakajima test setup is based on ISO12004 [25] with the modification that the test is conducted in a hydraulic press with a ram speed of 25 mm/s to reduce friction between punch and blank. The tool is the same one used for the Bulge test in the previous sub section. For Nakajima testing the punch is swapped to a spherical one in steel with diameter 100 mm. Strain data at the top surface of the sheet is evaluated and used for inverse modeling in this paper from:

  1. 1.

    FLC test with a 200 mm circular blank.

  2. 2.

    FLC test with a 125 mm wide dog bone sample.

  3. 3.

    LDH test with a 200 mm circular blank.

The FLC tests are carried out with lubrication conditions conforming to ISO12004 [25] to ensure a friction coefficient close to zero and the highest strains occurring close to the pole of the sample.

The LDH test is carried out with a cleaned and dry blank. This will yield an unknown friction level in the test. This test will be used to show that ambiguities occur in yield surface calibration if the friction is not known in the Nakajima test.

FE-models and discussion around UMAT verification

Two FE software are used for the simulations presented in this paper: LS-Dyna R9.3.0 and AutoForm R8.01. LS-Dyna is used for the implementation of BBC05 and all of the inverse modeling. AutoForm is only used to verify the implementation of BBC05 in LS-Dyna.

LS-Dyna

All the simulations of the Nakajima tests are carried out with the explicit FE code LS-Dyna. Double symmetry is utilized in all the models. Die velocity is 500 mm/s with a time step of 1.2 ∗ 10−6s and a resulting mass scaling ratio of ≈4 in all the models. The blanks are meshed with 2 mm rectangular Belytschko-Tsay elements with 5 integration points across the thickness. No mesh adaptivity occurs in the evaluated area of the blank. The tools are meshed with rigid body meshes. The contacts between blank and tool surfaces are of the type FORMING_ONE_WAY_SURFACE_TO_SURFACE.

AutoForm

AutoForm is an implicit FE-code where BBC05 with non-integer exponent already is implemented. The software is used as a reference for the implementation carried out in LS-Dyna in this paper. The model can be executed with arbitrary ram velocity since there are no time dependencies present. Double symmetry is utilized for the same tool and blank geometries as in the LS-Dyna model. The blank is meshed with the AutoForm triangular EPS11 elements. No mesh adaptivity occurs in the evaluated area of the blank.

Model differences and verification of BBC05 UMAT implementation

The explicit LS-Dyna model and the implicit AutoForm model gives very similar results for the prediction of strain and punch force, thereby verifying the implementation of the BBC05 model (BBC05 and YLD2000 gives equivalent results in LS-Dyna). The major and minor strain levels are within ≈0.01 strain units from each other. There is always differences when solving the same equation in different FE solvers. Here, the coulomb friction level needs to be adjusted 0.01 − 0.02 units when comparing AutoForm and LS-Dyna models. After localization, the LS-Dyna model has strains increasing faster than the AutoForm model, and the force curves drops a couple of millimeters earlier in punch depth. This can have several causes, e.g. element behavior, effects from high velocity and mass scaling in the explicit model, differences in implementation in the two FE software, refinement level of the models, mesh adaptivity, et.c. However, this paper is only using results before localization of the material where the two models behaves very similar. Users must always consider these effects when transferring results and knowledge from one FE code to another, e.g. from LS-Dyna to AutoForm which is the production code at Volvo Cars. The purpose of using the LS-Dyna Nakajima model is to be able to connect the process of material characterization to software such as MATLAB, and tune the model for speed while retaining a reliable level of strain and force predictions. Model speed is crucial when running large number of simulations, such as in optimization routines. A MATLAB script executing the LS-Dyna solver was used for creating Fig. 9.

Calibration of yield surface

Nakajima tests are selected as input for tuning the exponent k. Strains are considered in two sections of the blank: one in the rolling direction and one perpendicular to the rolling direction (i.e. the transverse direction). The inverse modeling of the Nakajima tests also measures and incorporates the punch force.

Friction in Nakajima tests causing ambiguities in yield surface calibration

Friction is an unknown parameter in the Nakajima simulations and must be identified through inverse modeling or advanced friction modeling, unless it holds true that friction can be neglected when tests are carried out according to ISO12004-2 [25]. Furthermore, if the exponent k in the BBC05 model is not set to a fixed value depending on material type but instead allowed to vary, two parameters need to be identified through inverse modeling.

The LDH test is here used to demonstrate ambiguities that occur in the yield surface calibration if friction is an unknown parameter. The results in sub sections Calibration of BBC05 with non-integer exponent and Coulomb friction and Calibration of BBC05 with integer exponent and Coulomb friction demonstrates that it is possible to calibrate the yield surface using the LDH results with different settings for friction and k yet achieve the same level of accuracy in the strain and force prediction. Of course, only one setting is physically correct; however, it cannot be found simply by considering strain levels and punch force in the LDH test. The plot in Fig. 9 visualizes these challenges. The surface is the fit between simulations and experiments for a combination of major and minor strain in the LDH test. The root mean square (RMS) values for the curves support the calculation of the fit, and both major and minor strains are normalized with the maximum strain values from the experimental measurements. The figure illustrates a clear linear dependence between the exponent k and the Coulomb friction coefficient. The plot in Fig. 9 was generated in MATLAB based on simulations with LS Dyna.

The section that is presented in Fig. 9 and the following two sub sections reflects a point in the LDH test at which the maximum plastic strain in the sheet is approximately 0.36, which is slightly after the tensile test curve but still on the bulge curve. The punch depth is 34.11 mm, before necking occurs.

Calibration of BBC05 with non-integer exponent and coulomb friction

An accurate fit of strains and forces in the LDH test for BBC05 with non-integer exponent can be obtained with the following parameters (for a discussion about differences between LS-Dyna and AutoForm, see sub section Model differences and verification of BBC05 UMAT implementation):

Yield surface exponent: k = 2.6.

Coulomb friction coefficient: μ = 0.20 (0.22 in AutoForm)

Calibration of BBC05 with integer exponent and coulomb friction

An accurate fit of strains and forces in the LDH test for BBC05 with integer exponent can be obtained with the following parameters (for a discussion about differences between LS-Dyna and AutoForm, see sub section Model differences and verification of BBC05 UMAT implementation):

Yield surface exponent: k = 3

Coulomb friction coefficient: μ = 0.16 (0.18 in AutoForm)

Separation of friction level and yield surface exponent by FLC simulations

Identifying k by simulations of the LDH test is clearly inhibited by the unknown parameter levels. Additional simulations and experiments is needed to fully separate the friction level, they should preferably exhibit different strain paths and have known friction levels. Hypothetically, this is already available in the FLC tests that are carried out for material characterization at Volvo Cars. The LDH test is carried out with friction, but the FLC tests should, according to ISO12004-2 [25], be prepared so that the test is nearly frictionless in the evaluated area. This should enable friction to be neglected in the inverse modelling process thereby making k the only unknown variable.

Several different methods to achieve a nearly frictionless test have been suggested in the literature, often consisting of one or several sheets of teflon, polyurethane, and additional lubrication between sheet and blank [25,26,27]. These sheets can be several millimeters thick. It is therefore important to measure the diameter on the top of the dome during the drawing of the blank. This can be done in the ARAMIS software by fitting a sphere to the top of the blank. The obtained diameter should be reduced with the blank thickness, and this new diameter is the punch diameter that shall be used in the simulation. It should also be checked that this diameter is relatively constant throughout the interval in the experiment that is used in the inverse modeling.(Figs. 5,6,7,8,9,10,11,12,13,14,15 and 16)

Fig. 5
figure 5

Final plastic hardening curve

Full size image
Fig. 6
figure 6

A) Nakajima/Bulge Die. B) Camera View. C) Example of results in Aramis Software

Full size image
Fig. 7
figure 7

LS-Dyna model during forming operation

Full size image
Fig. 8
figure 8

Example of a Nakajima simulation together with sections in rolling direction (RD) and Transverse Direction (TD) used for evaluation of strains

Full size image
Fig. 9
figure 9

Fit between experimental and simulated strains for LDH test (lower value = better fit). The red markers corresponds to settings in sub section Calibration of BBC05 with non-integer exponent and Coulomb friction and Calibration of BBC05 with integer exponent and Coulomb friction

Full size image
Fig. 10
figure 10

Measured and simulated strains at section in rolling direction through center of blank

Full size image
Fig. 11
figure 11

Measured and simulated strains at section in transverse direction through center of blank

Full size image
Fig. 12
figure 12

Measured and simulated punch force in LDH test

Full size image
Fig. 13
figure 13

Measured and simulated strains at section in rolling direction through center of blank

Full size image
Fig. 14
figure 14

Measured and simulated strains at section in transverse direction through center of blank

Full size image
Fig. 15
figure 15

Measured and simulated punch force in LDH test

Full size image
Fig. 16
figure 16

125 mm FLC Sample utilizing double symmetry

Full size image
Fig. 17
figure 17

Punch friction study of FLC 125 mm in Rolling Direction

Full size image
Fig. 18
figure 18

Punch friction study of FLC 125 mm in Transverse Direction

Full size image
Fig. 19
figure 19

Binder and die friction study of FLC 125 mm in Rolling Direction

Full size image
Fig. 20
figure 20

Binder and die friction study of FLC 125 mm in Transverse Direction

Full size image
Fig. 21
figure 21

Exponent k study of FLC 125 mm in Rolling Direction

Full size image
Fig. 22
figure 22

Exponent k study of FLC 125 mm in Transverse Direction

Full size image
Fig. 23
figure 23

Exponent k study of FLC 200 mm in Rolling Direction

Full size image
Fig. 24
figure 24

Exponent k study of FLC 200 mm in Rolling Direction

Full size image
Fig. 25
figure 25

Final FLC 125 mm in Rolling Direction

Full size image
Fig. 26
figure 26

Final FLC 125 mm in Transverse Direction

Full size image
Fig. 27
figure 27

Final FLC 200 mm in Rolling Direction

Full size image
Fig. 28
figure 28

Final FLC 200 mm in Rolling Direction

Full size image

The FLC samples are prepared from circular cutouts from the blanks according to ISO12004-2. One FLC sample is a circular sample of 200 mm corresponding to equi-biaxial tension. The other is a dog bone geometry according to [25] with the width 125 mm across the rolling direction, this sample corresponds to plane strain in the transverse direction of the blank. The samples have different widths to create different strain paths in the tests. Three simulation studies are carried out and compared to experimental results. Double symmetry is utilized in the FLC simulations. The FLC tests are evaluated at a punch depth at plastic strain close to the end of the hardening curve from the tensile test.

FLC 125 mm: Evaluated punch depth: 29.34 mm.

FLC 200 mm: Evaluated punch depth: 25.18 mm.

Friction in FLC simulations

First, two studies are carried out for the 125 mm FLC sample with a constant k value of 3. The first study has a constant friction on binder and die, the second one has a constant punch friction. These studies verifies the assumption that the friction level is 0.00 on the punch in the FLC tests carried out according to ISO12004-2. The major strain is highly dependent on the punch friction level, the only way to match experimental values is to set the Coulomb friction to 0.00. The binder and die friction mostly influences the minor strain levels, if it drops too low it also has an influence on the major strains. This is visualized in figs. 17,18,19,20,21, and 22.

Conclusions for friction levels are:

  • Punch friction: 0.00

  • Binder and Die friction: 0.06–0.07

There is a small tradeoff between major and minor strain prediction if 0.06 or 0.07 are used as binder and die friction. 0.06 are used for the remainder of the FLC simulations.

Exponent k in FLC simulations

A third study is carried out for the FLC 125 mm sample. It consists of a series of k from 2.5–3.5 where the punch friction is kept at 0.00.

Previously it has been concluded that the punch friction level in the FLC tests are very close to 0.00. k and punch friction are also nearly independent of each other in the studied range. k is influencing the strains on top of the punch, and punch friction influences the level of major strain. There is a small error in minor strain levels but it does not come from a wrong friction level on the punch since the minor strain is independent of the punch friction. Volvo Cars considers that a reliable calibration of a material gives strain levels within 0.01 units from the experimental values. The minor strains are within this interval.

After the three studies on the 125 mm sample is completed, another study is performed on the 200 mm FLC sample. The physical lubrication systems are identical in the 125 mm test and the 200 mm test. The punch friction in the 200 mm FLC simulation is therefore 0.00, and the binder and die friction is 0.06. Thereby, the only free variable in the 200 mm FLC simulations is k. Simulations with 3 different k are carried out and presented in Figs. 23 and 24.

Results: Final FLC simulations

Following conclusions can been made from the FLC simulations:

  • FLC punch friction level: 0.00

  • FLC binder and die friction level: 0.06–0.07

  • k is now the only free variable in the FLC simulations

The best fit of strains in the FLC samples is achieved with k ≈ 2.5. This is depicted in Figs. 25, 26, 27 and 28.

Conclusion and future work

The implementation of the BBC05 model with non-integer exponent k produces a flexible yield criterion. It is described how to implement BBC05 in FE software as a user subroutine, and calibration demonstrate higher accuracy and flexibility in strain predictions when the formulation with non-integer exponent is employed. However, when the BBC05 model with non-integer exponent was used for the calibration of the yield surface, uncertainties can emerge if friction is an unknown parameter. Previous research has reported the influence of friction in Nakajima, and similar tests [26], which warrants further investigation and consideration.

Two simulations of Nakajima samples are proposed for the inverse modeling process to separate the variables k and friction. The two Nakajima samples corresponds closely to the points of equi-biaxial tension and plane strain in transverse direction. The most reliable strain prediction is achieved with k=2.5, thereby also demonstrating the need for flexible yield criteria with non-integer exponent. It should also be noted that using LDH inverse modeling for BBC05 with integer exponent only yields a reliable prediction if the final exponent really is an integer number.

The authors suggest further exploration of friction in the Nakajima tests. Here, the experiments used for inverse modelling were carried out with friction conditions conforming to ISO12004-2 [25]. This should enable friction to be neglected in the simulations as noted in previous research, e.g. [19, 20]. The suggested process for inverse modelling confirms that the Coulomb friction coefficient shall be set to 0, the method should also detect and determine the friction punch level for cases where friction is unneglectable. Further studies should be carried out where different controlled combinations of lubrication, Teflon, and soft plastic sheets are used to control the friction. Preferably the amounts of liquid lubrication should be carefully measured and characterized, then used in the inverse modelling process together with FE representations of the Teflon and plastic sheets. An alternative route is to experimentally find the Coulomb friction in LDH tests, thereby yielding a test where k is the only unknown. This can however require friction tests for every material calibration performed.

A question that should be considered is that simulations predicts nearly equivalent major and minor strain at the pole of the Nakajima 200 mm sample, this is not the case in the DIC measurements where there seems to be some noise in strains on the pole. Can this be an effect of the soft and thick lubrication system described in ISO12004-2 [25]?

An indication that the friction actually is ≈0 is that the maximum strain occurs on the top of the dome in Nakajima tests, this is actually not the case for the 200 mm Nakajima samples presented in this paper. The cause of this should be investigated since the proposed method for inverse modeling confirms a friction level of ≈0. This behavior is also predicted in the simulations, so it could be an effect of the curved geometry and/or the shape of the yield surface.