Introduction

Dual phase (DP) steels are used in structural and safety parts of automotive vehicles such as longitudinal beams and reinforcements. These steels have enhanced strength and energy absorption capability thanks to their composite-like microstructure, which combines a hard martensite phase with a soft ferrite phase. The microstructure is obtained by specific control of the selection of alloying elements and the thermomechanical processes, leading to various grades with different microstructural and mechanical properties [1]. The microstructure of DP steels has important effects not only on the strength and ductility, but also on the anisotropic stress–strain behaviour observed under strain path changes [2–7].

Considering the complex load paths involved in the forming processes of automotive parts, numerical simulation of the stress–strain behaviour of DP steel under strain path changes is critical to determine the mechanical properties of the manufactured parts. However, the mechanical response of DP steel exhibits a complex dependence on the microstructure, where the chemical composition and microstructure characteristics of the ferrite and martensite phases are the main influencing parameters.

The effects of the alloying elements in solid solution depend on whether the element is interstitial or substitutional [8–17]. The interstitial elements such as carbon have stronger effects on the mechanical properties than substitutional elements due to their strong interaction with dislocations. It is often assumed that each element generates an increase in the yield strength proportional to its concentration and independent of the presence of other elements [8–17]. Due to the diffusion limit in the ferrite lattice, the carbon content in the martensite phase is much higher than the nominal content in the DP steel depending on the martensite volume fraction. An increased amount of carbon in a DP steel will therefore contribute to increasing the contrast in the stress–strain behaviour between the ferrite and martensite phases. Although, the effects of carbon on yield stress, hardness, and ultimate tensile strength of DP steels have been studied [4, 14–17], the effects of carbon on the contrast between the stress–strain response of the ferrite and martensite phases and on the residual stresses after deformation have received less attention.

The stress–strain response of DP steels is also affected by the volume fraction of the ferrite and martensite phases as well as the average size of ferrite grains and the distribution of martensite islands. An increased fraction of martensite increases the yield stress and ultimate tensile strength and most likely reduces the tensile elongation [4, 7, 15, 16, 18, 19]. The yield stress and ultimate tensile strength are improved by refining the ferrite and martensite microstructure [19–23]. The ferrite and martensite microstructure also affects the storage and patterning of dislocations and the internal stresses resulting from the contrast between the individual stress–strain response of the phases [4, 24]. Advanced experimental studies based on in situ testing methods, such as nanoindentation and nanopillars compression, have been conducted to assess the stress–strain response of the ferrite and martensite phases [25–27]. To the same end, in situ characterization, such as high-energy X-ray diffraction, has been combined with polycrystalline plasticity models [28, 29]. However, the measurement of the phase responses in DP steels is still a challenging task.

Experimental studies on the effects of the microstructure on the stress–strain response of DP steels under strain path changes have been limited to laboratory-enabled load paths [2–7, 24, 30]. Due to the limitations of the experimental test methods, there is a need for a systematic numerical approach which enables the prediction of microstructural effects under strain path changes. Simplified analytical approaches based on the law-of-mixture are not sufficiently accurate for this purpose [31–33]. Advanced analytical approaches, based for instance on the self-consistent model, have been used successfully to predict the response of DP steels under certain load paths [16, 34–37]. However, the Eshelby theory of an ellipsoidal inclusion embedded in an infinite medium has certain limitations on the representation of the microstructure and the phase interaction. To overcome this limitation, the finite element based representative volume element (RVE) approach has been used to account for the effects of the phase distribution on the mechanical properties of DP steels [23, 38–47]. In this approach, the DP steel microstructure is explicitly represented by the finite element mesh, which is divided into two sets of elements having the properties and spatial distributions of the ferrite and martensite phases.

An accurate representation of the contrast between the individual responses of ferrite and martensite phases is key to replicate the stress and strain partitioning in the microstructure by the RVE approach. Currently, the individual responses of ferrite and martensite might be determined by either empirical models or physically based models. While the empirical models fairly well predict the response of DP steel in specific conditions [5, 15, 16, 43, 44], they depend heavily on calibration tests relevant for the given application. The physically based models, which are formulated based on the dislocation theory, are less dependent on calibration tests. In this category, the crystal plasticity-based models show the greatest potential. These models determine the stress in each phase from the stresses in the grains considering single crystal plasticity by dislocations slip on active slip systems and crystal lattice rotation [28, 29, 36, 40]. Crystal plasticity-based simulations are CPU time intensive, which limits the physical size of the RVE and the finite element model. An alternative approach is to apply the simplified dislocation-based phase model [10–13] in the finite element based RVE model. The flow stress for each phase is computed as a function of the content of substitutional and interstitial solute elements, the grain size (i.e., the Hall–Petch effect) and the dislocation density, which is controlled by the storage and annihilation of dislocations. In these models, the variables controlling the storage and annihilation of dislocations in the ferrite phase are related to the average size of the ferrite grains [10–13, 38, 39, 41]. Due to the complex hierarchical microstructure of the martensite phase, the link between the dislocation mechanisms and characteristics of the martensite microstructure is not fully established. As a result, the model parameter controlling the storage of dislocations in the martensite phase is typically fitted to the experimentally obtained response of DP steels [38, 39, 41, 42, 46]. In this work, an extended dislocation-based phase model of martensite phase is proposed. The model parameter controlling the storage of dislocations is related to the carbon content in the martensite phase, which appears to be the most important parameter affecting the strength of martensite due to carbon’s strong interaction with dislocations. Particular attention is paid to the identification of model parameters and their effects on the contrast between ferrite and martensite responses, and further on the effects of pre-straining on the yield locus and the stress–strain behaviour of DP steels.

The paper is organized as follows. “Experiments” section introduces the material composition and microstructure characterization of DP and martensitic steels considered for the calibration of the phase-models and the establishment of RVEs for the investigated DP steels. “RVE-based modelling approach” section introduces the physically based constitutive model for the two phases, including the extended dislocation-based model of the martensite phase. The microstructure-based modelling framework of DP steels subjected to different loading paths is also presented in this section. In “Identification of material parameters” section, the calibration procedure for the extended phase model for martensite is presented. The effects of the alloying elements and the phase microstructure on the contrast between the ferrite and martensite responses are analysed. In “Numerical study” section, the modelling framework is applied to study the effects of pre-straining on the yield locus of DP steels and their stress–strain behaviour under reversed loading paths. A short summary of the results and perspectives for further development are given in “Conclusions” section.

Experiments

The steel materials studied in this work are four commercial grades of DP steels and three commercial grades of martensitic steels. The selected grades have different chemical compositions and microstructural characteristics. Tensile tests were performed on martensitic steel to calibrate the extended phase model of martensite, and on DP steels to evaluate the phase models for ferrite and martensite and the microstructure-based modelling framework.

Materials and microstructure

The chemical compositions of the DP steels and the martensitic steels are given in Tables 1 and 2, respectively. The DP steel grades are referred to as 500DP, 600DP, 800DP and 1000DP, whereas the martensitic steels are denoted 1200 M, 1400 M and 1700 M. These materials were supplied by SSAB as cold rolled sheet materials developed for the automotive industry. The sheet thickness was 2 mm for 500DP and 1 mm for all the other materials.

Table 1 Volume fraction of martensite phase, average grain size of ferrite phase, weight fraction of alloying elements in the DP steels, and weight fraction of carbon in martensite phase computed by Eq. (1)
Full size table
Table 2 Weight fraction of alloying elements in the martensitic steels
Full size table

The average grain size of the ferrite phase and the volume fraction of the martensite phase representative of commercial DP steel grades, which are similar to those investigated here, are reported in [33]. An exception is an additional amount of Nb for 600DP, which might result in a slightly less refined microstructure. However, the potential change in the microstructure was assumed to be negligible, and the microstructure measurements reported in [33] were used for all the investigated DP steels. Figure 1 shows scanning electron microscope (SEM) images of the microstructures of these DP steels grades from [33]. It should be noted that the microstructure of the investigated DP steels might deviate from the microstructure presented in Fig. 1. However, the effect of the potential deviation in microstructure for each DP steel was assumed negligible compared to the difference in microstructure between the four DP steel grades.

Fig. 1
figure 1

Scanning electron microscope (SEM) images representative of the microstructure of the investigated DP steels from [33]. The martensite phase appears as light areas and the ferrite phase as dark grey areas

Full size image

The substitutional elements were assumed to be equally distributed in the ferrite and martensite phases with contents equal to the nominal content in DP steels [10, 11, 39, 41]. The carbon content in the ferrite phase was approximated to 0.001%, which is the solubility limit of carbon in ferrite at room temperature [48]. The carbon content in the martensite phase was computed by the law of mixture from the nominal content of carbon:

$${C}_{C}^{DP}={V}_{m}{C}_{C}^{m}+\left(1-{V}_{m}\right){C}_{C}^{f}$$
(1)

where \({C}_{C}^{DP}\), \({C}_{C}^{f}\) and \({C}_{C}^{m}\) are the weight fractions of carbon in the DP steel, the ferrite phase, and the martensite phase, respectively, and \({V}_{m}\) is the volume fraction of the martensite phase. Table 1 shows the weight fractions of carbon in the martensite phase computed by Eq. (1) for each of the DP steels.

Tensile tests

The tensile tests were performed on flat specimens prepared from the DP and martensitic steel sheets by machining. The specimen geometry is illustrated in Fig. 2. Four specimens aligned with the rolling direction were prepared for each material. The tensile tests were performed at room temperature and quasi-static strain rate on a Zwick tensile machine with a 30 kN loadcell. Digital image correlation (DIC) was used to measure the displacement field on the specimen surface from which the tensile strains were calculated. A reference gauge length of 20 mm was used to determine the elongation and engineering strains.

Fig. 2
figure 2

Geometry of the flat tensile test specimen

Full size image

The engineering and true stress–strain curves for the DP steels and the martensitic steels are presented in Figs. 3 and 4, respectively. With one exception, the stress–strain curves show good repeatability until the ultimate tensile strength is reached. The exception is the set of stress–strain curves for 1700 M where scatter is observed before necking sets in. The average values of the yield stress measured at plastic strain \({\varepsilon }^{p}=0.5\mathrm{\%}\) (\({\mathrm{YS}}_{0.5}\)) and the ultimate tensile strength (UT) for each material are listed in Table 3 with the deviation computed by \(s/\sqrt{n}\), where \(s\) is the standard deviation and \(n\) is the number of tests.

Fig. 3
figure 3

Tensile tests on DP steels: (a) engineering stress–strain curves and (b) true stress–strain curves

Full size image
Fig. 4
figure 4

Tensile tests on martensitic steels: (a) engineering stress–strain curves and (b) true stress–strain curves

Full size image
Table 3 Measured yield stress (YS) and ultimate tensile strength (UT) of the DP steels and the martensitic steels
Full size table

RVE-based modelling approach

Single-phase model of ferrite and martensite

The phase model of ferrite and martensite is based on the single-phase steel model developed by Rodriguez et al. [10–13]. The ferrite and martensite phases are assumed to be homogeneous and isotropic materials and are modelled as von Mises materials with flow stress \({\sigma }_{f}\) as described in the following. The anisotropic response of DP steels under strain path changes is attributed to the residual stresses resulting from the contrast between ferrite and martensite responses. It should be noted that the ferrite and martensite phases develop substructural organizations of dislocations during deformation, leading to intra-phase kinematic hardening. In the current model, the effects of intra-phase dislocation structures were neglected, and the only effect of the contrast between the ferrite and martensite phases was accounted for.

The flow stress \({\sigma }_{f}\) for each phase is determined by:

$${\sigma }_{f}({\varepsilon }_{p})={{\sigma }_{y}+{\sigma }_{h}({\varepsilon }_{p})=\sigma }_{0}+{\sigma }_{s}+{\sigma }_{c}+{\sigma }_{d}+{\sigma }_{h}({\varepsilon }_{p})$$
(2)

where \({\sigma }_{y}\) and \({\sigma }_{h}\) represent the yield stress and the work hardening, respectively, and \({\varepsilon }_{p}\) is the equivalent plastic strain. The yield stress \({\sigma }_{y}\) has contributions from lattice friction, \({\sigma }_{0}\), strengthening from substitutional solute elements, \({\sigma }_{s}\), strengthening from carbon, \({\sigma }_{c}\), and strengthening from microstructure refinement, i.e., the Hall–Petch effect, \({\sigma }_{d}\). The contribution from lattice friction \({\sigma }_{0}\)=77 MPa was considered for both ferrite and martensite phase.

Considering that each alloying element generates an increase in the yield strength proportional to its concentration and independent of the presence of the other elements [9, 10], the strengthening by substitutional solute atoms is described by:

$${\sigma }_{s}=\sum_{i=1}^{{n}_{e}}{k}_{i}{C}_{i}$$
(3)

where \({C}_{i}\) is the concentration of a specific alloying element \(i\) (in \(\mathrm{wt}. \%\)), \({k}_{i}\) is the corresponding proportionality constant, and \({n}_{e}\) is the number of substitutional elements. The proportionality constants from [10, 49] were used for both ferrite and martensite phases: \({k}_{Mn}=80 \mathrm{MPa}/\mathrm{wt}. \%\), \({k}_{P}=750 \mathrm{MPa}/\mathrm{wt}. \%\), \({k}_{Si}=60 \mathrm{MPa}/\mathrm{wt}. \%\), \({k}_{Cu}=80 \mathrm{MPa}/\mathrm{wt}. \%\), \({k}_{Ni}=45 \mathrm{MPa}/\mathrm{wt}. \%\) and \({k}_{Mo}=11 \mathrm{MPa}/\mathrm{wt}. \%\), where the subindex indicates the alloying element. The strengthening by interstitial elements is limited to the effect of the carbon content:

$${\sigma }_{C}={a}_{C}{C}_{C}+{b}_{C}$$
(4)

where \({C}_{C}\) is the concentration of carbon (in \(\mathrm{wt}. \%\)), i.e., equal to \({C}_{C}^{f}\) for the ferrite phase and \({C}_{C}^{m}\) for martensite the phase. For the ferrite phase, \({a}_{C}^{f}=5000 \mathrm{MPa}/\mathrm{wt}. \%\) and \({b}_{C}^{f}=0\) were used, as proposed by Rodriguez and Gutierrez [10, 12, 49]. For the martensite phase, Rodriguez and Gutierrez [10, 12, 49] proposed \({a}_{C}^{m}=3065 \mathrm{MPa}/\mathrm{wt}. \%\) and \({b}_{C}^{m}=-161 \mathrm{MPa}\). However, in this study the values of \({a}_{C}^{m}\) and \({b}_{C}^{m}\) were calibrated based on experimental tests on martensitic steels (see “Identification of material parameters” section).

As already mentioned, the Hall–Petch effect represents the strengthening from grain boundaries, i.e., the strength of the phase increases with decreasing grain size, and is given by

$${\sigma }_{d}=\frac{{K}_{HP}}{\sqrt{d}}$$
(5)

where \(d\) is the grain size and \({K}_{HP}\) is a constant. In the current study, we will only account for the Hall–Petch effect in the ferrite phase, and thus \(d={d}^{f}\) represents the size of the ferrite grains. The martensite in DP steels has a complex microstructure composed of four hierarchical levels, including prior austenite grains, packets, blocks, and laths. It is unclear which of these morphological characteristics that represent the size effect in the Hall–Petch relation. Therefore, the strengthening due to the Hall–Petch effect was omitted for the martensite phase.

The work hardening is modelled by the Taylor equation:

$${\sigma }_{h}=\alpha M\mu b\sqrt{\rho }$$
(6)

where \(\rho\) is the dislocation density, \(M\) is the Taylor factor, \(\mu\) is the shear modulus, \(b\) is the length of the Burgers vector, and \(\alpha\) is a numerical constant. The evolution of the dislocation density is governed by the Kocks-Mecking equation [50]:

$$\frac{d\rho }{d{\varepsilon }^{p}}=M\left(\frac{1}{bL}-k\rho \right)$$
(7)

where \(k\) is the dislocation annihilation constant and \(L\) is the dislocation mean free path. Integration of Eq. (7) is performed by assuming that dislocation–dislocation interaction can be neglected [10], leading to a strain-independent dislocation mean free path, and the dislocation density is given by

$$\rho=\frac1{bLk}\left[1-\exp\left(-Mk\varepsilon^p\right)\right]+\rho_0\exp\left(-Mk\varepsilon^p\right)$$
(8)

where \({\rho }_{0}\) is the initial dislocation density. It is here assumed that \({\rho }_{0}=0\), and then the work hardening is obtained from Eqns. (6) and (8) as:

$${\sigma }_{h}=\alpha M\mu \sqrt{b} \sqrt{\frac{1-\exp(-Mk{\varepsilon }^{p})}{kL}}$$
(9)

The hardening parameters of the ferrite phase are related to the average size of the ferrite grains according to [11]. For ultra-low carbon steels up to a grain size of about 10 µm, the effective dislocation mean free path is equal to the ferrite grain size [11], i.e., \({L}^{f}={d}^{f}\). The annihilation parameter is assumed to be inversely proportional to \({d}^{f}\) [10, 11]:

$${k}^{f}=\frac{{c}_{k}^{f}}{{d}^{f}}$$
(10)

where \({c}_{k}^{f}\) is a constant determined by fitting and can be approximated to 8 µm for low carbon steels with ferrite grain size \(5\;\mu m<{d}^{f}<20\;\mu m\) [12].

For the martensite phase, an increased amount of carbon is assumed to affect the dislocation storage in low-carbon lath martensite [51]. Due to carbon’s strong interaction with dislocations, the dislocation mean free path \({L}^{m}\) is assumed to be inversely proportional to the carbon content in the phase:

$${L}^{m}=\frac{{c}_{L}^{m}}{{C}_{C}^{m}}$$
(11)

where \({c}_{L}^{m}\) is the proportionality constant for lath martensitic steel, which is determined by fitting to experimental data and \({C}_{C}^{m}\) is wt.% of carbon in the martensite phase. This assumption is consistent with microstructural observations of the effects of carbon solute on dislocation motion in lath martensite. Niino et al. [51] reported that the carbon solute reduces the mobility of dislocations and increases their multiplication during deformation. The annihilation parameter of martensite \({k}^{m}\) is assumed to be constant independent on the carbon content.

Representative volume elements (RVEs) of DP steels

The RVE approach is employed to account for the composite-like microstructure of DP steels. The RVE was designed to include sampling of the microstructural heterogeneities generated based on statistical parameters representative of the phase distribution, i.e., the volume fraction of martensite and the average size of the martensite islands [41, 43]. Compared to the RVE generation processes based on direct translation of microstructural images [39, 40, 44, 45, 47, 52], the statistical generation process has less limitation on the size of the RVE (i.e., it is not limited to the size of the microstructure images) and the procedures for generating the three-dimensional RVE are simplified.

The statistical generation involves a three-step partitioning procedure.

  • In the first step, a cuboid part with geometry of \(0.1\times 0.1\times 0.05 {\mathrm{mm}}^{3}\) is discretized into a finite element mesh using fully integrated cubic elements. The total number of the elements in each RVE is 62500. Each element has a unique reference associated to a set of parameters including the connectivity of its nodes and the coordinates of its centre.

  • In the second step, the elements are partitioned into arbitrary domains representing a periodic aggregate of equiaxed grains. The partitioning process is based on Voronoi tessellations and starts with the generation of nucleation germs representing the centre of the equiaxed grains [53]. The number of the nucleation germs was determined by \(N=V/{d}_{m}^{3}\), where \(V\) is the volume of the cuboid part and \({d}_{m}\) is the average size of martensite islands. Each element in the finite element mesh was assigned to the closest germ, where the distance between the centre of the element and the nucleation germ is determined by:

    $${d}^{en}=\sqrt{\sum_{i=1}^{3}{\left(\mathrm{min}\left(\left|{x}_{i}^{e}-{x}_{i}^{n}\right|,\left|{l}_{i}^{RVE}-\left|{x}_{i}^{e}-{x}_{i}^{n}\right|\right| \right)\right)}^{2}}$$
    (12)

    where, \({x}_{i}^{e}\), \({x}_{i}^{n}\) and \({l}_{i}^{RVE}\) are respectively the coordinates of the element centre, the coordinates of the nucleation germ, and the length of the RVE in the \({x}_{1}\)-, \({x}_{2}\)- and \({x}_{3}\)-directions. The lengths of the RVE are introduced in the second term of the minimum function to ensure the periodicity of the generated microstructure.

  • In the third step, the grains were randomly partitioned into two sets representing the martensite and ferrite phases. The partitioning process is controlled by the volume fraction of the martensite. The elements in each set are assigned to the associated phase reference. The number of elements assigned to each phase is adjusted when necessary to reach the exact value of the volume fraction. This adjustment did not affect the average size of martensite islands nor their distribution.

Figure 5 shows the RVEs of the investigated DP steels. In these RVEs, the microstructure periodicity is ensured by the periodicity of the grain aggregate generated from the second step. As can be observed, the generation process leads to step-like boundaries between the ferrite and martensite phases. Although these boundaries might affect the local response of the interfaces compared to smooth boundaries, their effects on the global response of the RVE are deemed negligible when a sufficiently refined mesh is applied.

Fig. 5
figure 5

RVEs of the four DP steels

Full size image

Finite element modelling of load paths

The overall or average stress state in the RVE is obtained by the finite element method as a volume average of the stresses over the individual phases. Each phase is explicitly discretized into \({n}_{e}^{p}\) elements with \({n}_{i}^{p}\) integration points. The Cauchy stress tensor \({{\varvec{\sigma}}}_{i}^{p}\) is computed for each integration point based on the single-phase model, and then the Cauchy stress tensor \({{\varvec{\sigma}}}^{p}\) for each phase is determined by:

$${{\varvec{\sigma}}}^{p}=\sum_{e=1}^{{n}_{e}^{p}}{f}_{e}^{p}\sum_{i=1}^{ {n}_{i}^{p}}{f}_{i}^{p}{{\varvec{\sigma}}}_{i}^{p}$$
(13)

where \({f}_{i}^{p}\) is the volume fraction of integration point \(i\) of element \(e\) and \({f}_{e}^{p}\) is the volume fraction of element \(e\) in phase \(p\). The average stress state of the RVE is computed by:

$${\varvec{\sigma}}=\sum_{p=1}^{2}{f}^{p}{{\varvec{\sigma}}}^{p}$$
(14)

where \({f}^{p}\) is the volume fraction of phase \(p\).

The load paths in the finite element simulations are created by applying displacement boundary conditions on the RVE as illustrated in Fig. 6 [54, 55]. The RVE has four master nodes located at the corners and six reference surfaces representing its external surfaces. Master node 1 is at the origin of the coordinate system applied for the RVE and is fixed in space, i.e., all degrees of freedom are set to zero. A specific load path is obtained by controlling the displacements of master nodes 2, 3 and 4, see Fig. 6a. The periodicity of the boundary conditions applied on the RVE is ensured by controlling the displacement of the nodes on the reference surfaces. The displacement vectors of each couple of symmetric nodes located on two opposite reference surfaces are related by:

$${{\varvec{u}}}^{b}-{{\varvec{u}}}^{a}={{\varvec{u}}}^{R}-{{\varvec{u}}}^{1}$$
(15)

where \({{\varvec{u}}}^{a}\) and \({{\varvec{u}}}^{b}\) are the displacement vectors of two symmetric nodes, such as nodes \(a\) and \(b\) in Fig. 6a, located on two opposite reference surfaces, \({{\varvec{u}}}^{R}\) is the displacement vector of the master node located in the same reference surface as node \(b\), and \({{\varvec{u}}}^{1}=\bf0\) is the displacement vector of the master node located at the origin. Equation (15) relates the displacements of the nodes on each pair of opposite faces to the displacement of the corresponding reference nodes, leading to periodic deformation of the RVE faces as demonstrated in [55].

Fig. 6
figure 6

Control of nodal displacements a) to create periodic boundary conditions, see Eq. (15), b) to compute the yield locus, c) to simulate forward-reverse shear loading, and d) to simulate tension–compression loading

Full size image

The RVE was used to calculate the yield loci of the DP steels, either with or without pre-straining in uniaxial tension, and to study the stress–strain behaviour of the DP steels under uniaxial tension–compression loading and forward-reverse shear loading. To calculate the yield locus, in-plane stress states were first generated according to the method established by Saai et al. [54, 56]. Only stress states at yielding with \({\sigma }_{12}=0\) were computed by the RVE model and plotted as discrete points on the yield locus [43]. Each yield locus was determined by four sequences of simulations, in which the nodal velocities of master node 2 and master node 3 were controlled:

$$\begin{array}{l}\begin{array}{c}\mathrm{Sequence}\;1:v_1^2=v_0\;\mathrm{and}\;v_2^3={\rho v}_0\\\mathrm{Sequence}\;2:v_1^2={-v}_0\;\mathrm{and}\;v_2^3={\rho v}_0\end{array}\\\mathrm{Sequence}\;3:v_2^3=v_0\;\mathrm{and}\;v_1^2={\rho v}_0\\\mathrm{Sequence}\;4:v_2^3={-v}_0\;\mathrm{and}\;v_1^2={\rho v}_0\end{array}$$
(16)

where \({v}_{i}^{2}\) and \({v}_{i}^{3}\) are the nodal velocities of master nodes 2 and 3 in the \({x}_{i}\)-direction, \({v}_{0}\) is a constant reference velocity, and \(\rho \in [-1, -0.6,-0.2, 0.2, 0.6, 1]\) is a proportionality factor that determines the load path. The nodal velocity components \({v}_{1}^{2}\) and \({v}_{2}^{3}\) are illustrated in Fig. 6b. Each simulation within a sequence gives a stress path and an associated point on the yield locus. Each sequence gives points at yielding that lie within one of the quadrants of the yield locus. It was carefully checked that the out-of-plane stress components \(({\sigma }_{33},{\sigma }_{23},{\sigma }_{31})\) computed with these boundary conditions were negligibly small compared to the in-plane stress components \(({\sigma }_{11},{\sigma }_{22},{\sigma }_{12})\). The yield locus was computed for each DP steel either with or without pre-straining in uniaxial tension. In the case without pre-straining, the RVE represents the initial microstructure of the material, while in the case with pre-straining, the RVE is first subjected to uniaxial tension in the \({x}_{1}\)-direction up to a total strain of \(0.0315\), and then unloaded before the computation of the yield locus.

The stress states at yielding were determined for each RVE at equivalent amount of plastic work per unit volume. The reference plastic work per unit volume is denoted \({W}_{p}^{0.x\%}\) which defines the plastic work per unit volume at a plastic strain of \(0.x \%\) in uniaxial tension along the reference direction. The reference direction was taken to be the \({x}_{1}\)-direction in Fig. 6a. For each RVE without pre-straining, the stress states at yielding were computed at two reference values of plastic work per unit volume, \({W}_{p}^{0.2\mathrm{\%}}\) and \({W}_{p}^{0.5\mathrm{\%}}\). For the pre-strained RVEs, the stress states at yielding were computed at an incremental plastic work per unit volume after pre-straining \(\Delta {W}_{p}\) equal to \({W}_{p}^{0.x\mathrm{\%}}\), where \(0.x\mathrm{ \%}=0.2\mathrm{\%} \; \mathrm{ or }\; 0.5\mathrm{\%}\) is the reference plastic strain in uniaxial tension along \({x}_{1}\)-direction after pre-straining. To investigate strain path change, forward-reverse shear and cyclic uniaxial tension–compression simulations were performed with each RVE. The boundary conditions used to simulate forward-reverse shear loading and tension–compression loading are illustrated in Fig. 6c and d, respectively.

The FE analysis was carried out using ABAQUS Standard (Dassault Systemes SIMULIA Corp). Python scripts were developed to generate the RVE mesh including nodes and elements sets, write periodicity equations and boundary conditions, and manage the simulation sequences for the computation of the yield locus. Python scripts were also developed for post-processing the output files and extract the yield locus and stress–strain response under strain path changes.

Identification of material parameters

The substitutional alloying elements were considered to have a uniform distribution in the ferrite and martensite phases with content equal to the nominal content in the DP steels. The effect of these elements on the yield stress is thus the same for the ferrite and martensite phases. The contribution from the substitutional elements to the yield stress \({\sigma }_{s}\) is calculated from Eq. (3) using the data provided in Table 1.

The tensile stress–strain curves of martensitic steels were used to calibrate the martensite phase model. This calibration process does not account for the difference between the microstructure of the martensite phase and the microstructure of martensitic steel. However, it is challenging to determine the stress–strain response of the martensite phase by experimental tests, as reported in the literature review [25, 26]. In the current model, the carbon content was assumed to be the key parameter influencing the strength of the martensite phase [14, 15, 18], and therefore its effects were calibrated based on the stress–strain response of martensitic steels with different carbon content.

For the ferrite phase, the contribution of the carbon content to the yield stress was determined by Eq. (4), adopting the constants \({a}_{C}^{f}\) and \({b}_{C}^{f}\) given by Rodriguez et al. [12, 49] and the solubility limit of carbon, i.e., \({C}_{C}^{f}=\) 0.001 wt.%. This contribution does not account for the effect of carbide particles, which are observed in the ferrite phase. The carbide particles will affect the calibration of the solute strengthening of the martensite phase, as will be discussed further in this section. For the martensite phase, the constants \({a}_{C}^{m}\) and \({b}_{C}^{m}\) in Eq. (4) were determined by fitting based on the tensile yield stresses of the three martensitic steels. Figure 7 shows the yield strengthening by carbon (\({\sigma }_{C}^{m}\)) and the ultimate tensile strength (UT) of the martensitic steels as function of the carbon content. The linear trendline in Fig. 7a gives \({a}_{C}^{m}=2442\mathrm{ MPa}/\mathrm{wt}.\mathrm{\%}\) and \({b}_{C}^{m}=417 \mathrm{MPa}\). From the figure it is seen that both the yield strengthening by carbon and the ultimate tensile strength tend to increase linearly with the carbon content.

Fig. 7
figure 7

a) Yield strengthening by carbon (\({\sigma }_{C}^{m}\)) and b) ultimate tensile stress (UT) of the three martensitic steels as a function of the carbon content

Full size image

The value of \({K}_{HP}^{f}\), which determines the strength of the Hall–Petch effect for the ferrite phase, has been found to be dependent on the alloying elements, where \({K}_{HP}^{f}\) is in the range of \(5 \mathrm{MPa}\sqrt{\mathrm{mm}}\) for interstitial free steels and 18 MPa \(\sqrt{\mathrm{mm}}\) for low carbon steels [57]. In the current work, \({K}_{HP}^{f}=14 \mathrm{MPa }\sqrt{\mathrm{mm}}\) was found to give the best fit to experimental responses of the investigated DP steels.

The hardening parameters of the ferrite phase, \({L}^{f}\) and \({k}^{f}\), were related to the average size of the grains, according to Gutierrez et al. [10, 11, 49]. For ultra-low carbon steels with grain size up to about \(10 \mathrm{\mu m}\), they reported that the effective dislocation mean free path can be replaced by the average grain size \({d}^{f}\). The value of \({k}^{f}\) was determined by fitting to experimental data, assuming \({k}^{f}\) to be inversely proportional to \({d}^{f}\), as expressed by Eq. (10). The proportionality constant \({c}_{k}^{f}\) \({=10}^{-5}\) m was applied [49]. These values of the hardening parameters have been applied in many numerical studies to simulate the stress–strain behaviour of the ferrite phase in DP steels [38, 39, 41]. In these studies, the hardening parameters of the martensite phase were calibrated to fit the stress–strain behaviour observed experimentally for the DP steels.

The hardening parameters of the martensite phase, \({c}_{L}^{m}\) and \({k}^{m}\), were calibrated based on the experimental measurements of the ultimate tensile stress (UT) of the martensitic steel. As illustrated in Fig. 7b, the ultimate tensile stress of the martensitic steels depends linearly on the carbon content. Based on this linearity, it was assumed that the dislocation mean free path is inversely proportional to the carbon content, cf. Equation (11), and that the annihilation parameter is independent of it. Fitting to the experimental data gave \({c}_{L}^{m}=0.0217 \mathrm{\mu m}\cdot \mathrm{wt}.\mathrm{\%}\) and \({k}^{m}=5.4\times {10}^{-5}\). Figure 8 shows the dislocation mean free path as predicted by Eq. (11) for the martensitic steels and the martensite phase in DP steels. Figure 9 compares the stress–strain response of martensitic steels as predicted by the phase model to the experimental measurements. The phase model predicts the response of martensitic steels with satisfactory accuracy, but the stress level is underestimated in the elastic–plastic transition.

Fig. 8
figure 8

Dislocation mean free path expressed by Eq. (11) as a function of carbon content for martensitic steels (crosses) and martensite phase in DP steels (dots)

Full size image
Fig. 9
figure 9

Equivalent stress versus equivalent plastic strain curves computed by the single-phase model for martensite compared to experimental stress–strain curves of martensitic steels

Full size image

The calibrated phase models were used in RVE simulations to calculate the uniaxial stress–strain curves of the DP steels. The results are shown in Fig. 10a. While the stress–strain behaviour is well predicted for 500DP, 600DP and 800DP, the deviation is substantial for 1000DP with the highest volume fraction of martensite. The most probable reason for the overestimation of the stress level for 1000DP is that the strength of the martensite is overestimated by the phase model, whereas the work hardening is well predicted. The carbon content in the martensite phase was estimated by a law-of-mixture, assuming that the carbon content in the ferrite was given by the solubility limit and neglecting the carbide particles, which have been observed in the ferrite phase [33]. The presence of carbides in the ferrite leads to an overestimation of the carbon content in the martensite and consequently the contrast between ferrite and martensite is exaggerated. Since the carbon content in the martensite phase was not measured, the constant \({b}_{C}^{m}\) in Eq. (4) was recalibrated based on the results in Fig. 10a, and the new value was found to be \({b}_{C}^{m}=117\mathrm{ MPa}\). The results of the RVE simulations with the recalibrated value of \({b}_{C}^{m}\) are shown in Fig. 10b, and it is clearly seen that the agreement between the simulation and experiment for 1000DP is improved, while for the other DP steels there are only minor changes in the simulation results.

Fig. 10
figure 10

Stress–strain curves of DP steels as predicted by the RVE simulations using a single-phase model of martensite with a) \({b}_{C}^{m}=417 \mathrm{MPa}\), as calibrated using the yield stress of the three martensitic steels, and b) \({b}_{C}^{m}=117 \mathrm{MPa}\) after adjustment to the measured stress–strain curve of 1000DP steel

Full size image

The equivalent stress versus plastic strain curves for the ferrite and martensite phases of the DP steels after the final calibration are shown in Fig. 11. While the ferrite phase exhibits work hardening that persists up to large strain, the work hardening of the martensite phase saturates at low strain. The difference in strength between the ferrite and martensite phases is mostly controlled by the higher carbon content of the martensite which increases the yield strength by solute strengthening.

Fig. 11
figure 11

Equivalent stress versus equivalent plastic strain curves of a) the ferrite phase and b) the martensite phase of the DP steels as predicted by the single-phase models after final calibration

Full size image

Numerical study

In-plane yield loci

The stress states at yielding were computed for 24 in-plane load paths giving the yield locus in the \({\sigma }_{11}-{\sigma }_{22}\) plane for \({\sigma }_{12}=0\). To evaluate the anisotropy, geometrically presented by the shape of the yield locus [54], the stress states at yielding were normalized by the uniaxial yield stress (\({\sigma }_{0}\)) along the \({x}_{1}\)-direction: \(\left|{\sigma }_{ij}\right|={\sigma }_{ij}/{\sigma }_{0}\). The normalized yield stresses were used to evaluate the effect of the microstructure of the DP steels and the pre-straining of the RVE. The yield loci of the DP steels computed with the RVEs without pre-straining are presented in Fig. 12a while the normalized yield loci are plotted in Fig. 12b. The yield loci are plotted for zero shear stress with yield points defined by \({W}_{p}\) equal to either \({W}_{p}^{0.2\mathrm{\%}}\) or \({W}_{p}^{0.5\mathrm{\%}}\). It appears that the RVE represents an isotropic material, and the shape of the yield locus fits accurately with the von Mises yield locus for all investigated DP steels. As reported in [4, 6, 58, 59], the anisotropic response of DP steel without pre-straining is mainly related to the crystallographic texture and the shape and non-homogenous distribution of the phases. In this work, the RVEs include only equiaxed grains (particles) described by an isotropic plasticity model, leading to an isotropic yield locus.

Fig. 12
figure 12

a) In-plane yield loci for zero shear stress and plastic work per unit volume \({W}_{p}\) equal to \({W}_{p}^{0.2\%}\) (left) and \({W}_{p}^{0.5\%}\) (right) at yielding, and b) corresponding yield loci normalized by the uniaxial yield stress along the \({x}_{1}\)-direction (\({\sigma }_{0}\))

Full size image

To introduce residual stresses and strains, the RVEs were pre-strained in uniaxial tension in the \({x}_{1}\)-direction up to a total strain of \(0.0315\). Then, they were relaxed to release the elastic deformations and to obtain a microstructure with residual stresses and strains. Figures 13 and 14 show the residual equivalent stresses and residual plastic strains resulting from the pre-straining, respectively. Figure 13 shows the effect of the microstructure on the gradient of stresses within each phase, which decreases with the increased volume fraction of martensite. The contribution of the martensite phase to the residual strain also depends on the microstructure, see Fig. 14. It is seen that 1000DP has the highest fraction of martensite elements exhibiting residual plastic deformation, while the 600DP and 800DP have higher maximum plastic strain localized in few elements on the phase interfaces. The fraction of martensite elements exhibiting residual plastic strain in 500DP, 600DP and 800DP increases with the volume fraction of martensite. However, the maximum values of stress and strain in the ferrite and martensite phases do not show a linear dependence on the volume fraction of the phase, see Figs. 13 and 14. This might be explained by the combined effects of the local response of the phase (Fig. 11) and the phase distribution.

Fig. 13
figure 13

Residual von Mises equivalent stress in the ferrite and martensite phases computed for each DP steel after pre-straining and unloading

Full size image
Fig. 14
figure 14

Residual von Mises equivalent plastic strain in the ferrite and martensite phases computed for each DP steel after pre-straining and unloading

Full size image

Figure 15 presents the yield loci of the DP steels computed with the RVEs after pre-straining in uniaxial tension along the \({x}_{1}\)-direction and up to a total strain of \(0.0315\). The yield loci are plotted for zero shear stress with yield points defined by \(\Delta {W}_{p}\) equal to either \({W}_{p}^{0.2\mathrm{\%}}\) or \({W}_{p}^{0.5\mathrm{\%}}\), where \(\Delta {W}_{p}\) is the incremental plastic work per unit volume of the RVE after the prestraining. It is seen that the yield loci are shifted along the \({\sigma }_{11}\)-axis because of the prestraining and that this shift is more pronounced for the yield loci defined by \(\Delta {W}_{p}={W}_{p}^{0.2\mathrm{\%}}\) than by \(\Delta {W}_{p}={W}_{p}^{0.5\mathrm{\%}}\).

Fig. 15
figure 15

In-plane yield loci for zero shear stress and incremental plastic work per unit volume after pre-straining \(\Delta {W}_{p}\) equal to \({W}_{p}^{0.2\%}\) (left) and \({W}_{p}^{0.5\%}\) (right) at yielding

Full size image

In Fig. 16, the yield loci presented in Fig. 15 were normalized by the yield stress \({\sigma }_{0}\) of the pre-strained RVEs in uniaxial tension along the \({x}_{1}\)-direction. They are compared to the normalized yield locus of the initial RVEs represented by von Mises yield locus, see Fig. 12b. As can be observed, the pre-straining of the RVEs introduces significant anisotropy. The yield loci of the pre-strained RVEs are shrunk and shifted compared to the yield locus of the RVEs without pre-straining. The shape of the yield loci and the amount of shrinking and shifting also vary depending on the DP steel microstructure. The 800DP and 600DP steels, with respectively 25% and 18% martensite, have higher anisotropy (due to the shrinking and shifting of the yield locus) than the 1000DP and 500DP steels, with respectively 50% and 15% martensite. The anisotropy results from the residual stresses in the RVE after pre-straining and depends on the reference value of the plastic work per unit volume at which the yield loci are computed, see Fig. 15. The RVE simulations predicted higher degree of shrinking of the yield loci defined by \(\Delta {W}_{p}={W}_{p}^{0.2\%}\) than of the yield loci at \(\Delta {W}_{p}={W}_{p}^{0.5\%}\), see Fig. 16.

Fig. 16
figure 16

Normalized in-plane yield loci for zero shear stress and incremental plastic work per unit volume after pre-straining \(\Delta {W}_{p}\) equal to \({W}_{p}^{0.2\%}\) (left) and \({W}_{p}^{0.5\%}\) (right) at yielding. The normalized yield locus of the RVE without pre-straining is plotted in grey, dashed line for comparison

Full size image

Stress–strain behaviour under load path change

The stress–strain behaviour of the DP steels under load path change was computed using the RVE approach with the single-phase models presented in “RVE-based modelling approach” section. Two cases were considered: forward-reverse shear loading and tension–compression cyclic loading. In the first case, each of the RVEs presented in Fig. 5 was deformed in forward simple shear to a total strain of 0.212, then followed by a reverse simple shear. Figure 17a shows the resulting stress–strain curves, where the shear stress magnitude is plotted as a function of the accumulated shear plastic strain. In the second case, each of the RVEs was subjected to a tension–compression cyclic load with a strain amplitude equal to 0.6%. The resulting cyclic stress–strain curves are plotted in Fig. 18a. For each case, the response of 600DP is compared to experimental data from Hérault et al. [5] on a similar DP steel in Figs. 17b and 18b, respectively. The numerical simulations predict a substantial Bauschinger effect for all DP steels in each case, which results from the contrast between the local response of the ferrite and martensite phases. To evaluate how the Bauschinger effect is influenced by the microstructure, the Bauschinger effect ratio (BER), in which the back stress is used and made dimensionless, is defined by [59,60,61]: \(\mathrm{BER}=\left(\left|{\sigma }_{RP}\right|-\left|{\sigma }_{YR}\right|\right)/\left|{\sigma }_{RP}\right|\), where \({\sigma }_{RP}\) is the stress at the reversal point and \({\sigma }_{YR}\) is the stress at yielding in reverse loading. The Bauschinger effect ratios of DP steels were computed for forward-reverse shear loading and tension–compression cyclic loading and are listed in Table 4. The RVE simulations demonstrate a clear dependency of the Bauschinger effect on the microstructure of the DP steels. The BERs increase with the martensite volume fraction, where 1000DP has the greatest Bauschinger effect. These observations are consistent with experimental observations by Watanabe et al. [60].

Fig. 17
figure 17

a) Stress–strain curves in forward-reverse shear predicted by RVE simulations for DP steels and b) predicted stress–strain curves for 600DP compared to experimental data from Hérault et al. [5]

Full size image
Fig. 18
figure 18

a) Stress–strain curves in cyclic uniaxial tension–compression predicted by RVE simulations for the DP steels and b) predicted curve of 600DP steel compared to experimental data from Hérault et al. [5]

Full size image
Table 4 Bauschinger Effect Ratio (BER) predicted by the RVE model for forward-reverse shear and tension–compression simulations performed on DP steels
Full size table

In Fig. 17b, the stress level after the elastic–plastic transition predicted by the RVE simulation for reverse shear loading reaches the stress level in simple shear loading. The experimental stress–strain curve of forward-reverse shear loading, however, shows a reduction in the stress level after load reversal compared to the stress–strain curve of simple shear. The RVE simulation also somewhat underestimates the Bauschinger effect for tension–compression cyclic loading compared to the experimental data (see Fig. 18b). The deviation might be explained by the effects of other microstructural features than the heterogeneity introduced by the contrast between ferrite and martensite responses, e.g., the evolution of dislocation arrangement within the phases [59, 62]. Gardey et al. [59] reported that these dislocation arrangements remained in most ferrite grains after an orthogonal strain path change, whereas they were partially dissolved under reverse loading. Notwithstanding, Figs. 17 and 18 show the important contribution of the contrast in strength between the ferrite and martensite phases on the elastoplastic transition after load reversal.

Conclusions

The yield stress and the strain hardening of the martensite phase in DP steel were related to the carbon content in the phase, which appears to be the most influencing parameter controlling the local response of the phase and its contrast with the local response of the ferrite phase. The contrasting behaviour between the ferrite and martensite phases and the microstructure characteristics of DP steel affect the residual stresses and residual strains in the RVE model leading to a complex anisotropic response under orthogonal and reverse loading.

The pre-straining of the RVEs introduces significant anisotropy represented by the shape of the computed yield loci. The yield loci of the pre-strained RVEs are shrunk and shifted compared to the yield loci of the RVEs without pre-straining. The amount of shrinking and shifting vary depending on the DP steel microstructure and the reference value of the plastic work per unit volume at which the yield loci are computed. The numerical simulations predict a substantial Bauschinger effect under tension–compression cyclic loading and forward-reverse shear loading. The contrast between the local response of the ferrite and martensite phases in the RVE model produces a significant Bauschinger effect which was found to increase with the amount of martensite phase. The single-phase models and the RVE model highlight the importance of the contrast between the local response of the ferrite and martensite phases on the structural response of DP steels under strain path changes.