**1. Introduction**

Dual-Phase (DP) steel, due to its higher energy absorption capacity and reduced weight, is used in the automotive industry to achieve simultaneous high strength and elongation goals [1–3]. Hard phase martensite laths embedded in the softer ferrite matrix are responsible for reinforcing the solid aggregate, while ductility is incorporated by the matrix [4]. Dual-phase steel is a suitable example of a multi-phase material because of the significant difference in the mechanical properties of its phases, and it has widespread usage in the automotive industry. This peculiar combination of hard and soft phases imparts desirable properties in the material, i.e., low 0.2% proof stress and a high work-hardening coefficient (n-value). Generally, during cold forming processes of drawing and stretching, higher n-values exhibit uniform global formability by avoiding local thinning [5]. However, practical aspects often reveal unexpected necking and local failure in complex forming during bending or flanging processes. Due to the heterogeneous microstructure of soft

**Citation:** Hussein, T.; Umar, M.; Qayyum, F.; Guk, S.; Prahl, U. Micromechanical Effect of Martensite Attributes on Forming Limits of Dual-Phase Steels Investigated by Crystal Plasticity-Based Numerical Simulations. *Crystals* **2022**, *12*, 155. https://doi.org/10.3390/ cryst12020155

Academic Editor: Wojciech Polkowski

Received: 20 December 2021 Accepted: 17 January 2022 Published: 21 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

ferrite and hard martensite, local straining causes unpredicted local abnormalities. The heterogeneity affects the micro-scale attributes of materials. Consequently, it influences the component scale's material properties, particularly the material damage behavior [6]. Therefore, it is of the utmost importance to investigate the relationship between the phases' heterogeneity and their microstructural attributes, especially martensite and ferrite fractions and their grain sizes [7].

Material formability is a fundamental mechanical property for vehicle bodies and structural members in the automotive sector. Therefore, accurate predictions of its indicators and reliable utilization are needed. For this purpose, FLDs are commonly used to estimate sheet metal formability. Its purpose is to predict safe, necking, and failure zones by adopting major and minor forming limit strain inside a diagram [8]. Experimentally, the Nakajima test detects these forming limit strains, which is used by performing a punch– die method on specimens with varying dimensions under different loading conditions. The test shows a higher accuracy for most materials, but it is expensive, slow, and relies on a complex specimen-shaping process. Therefore, FLDs are plotted using some wellestablished numerical models [9]. Tasan et al. [10] carried out nanoindentation experiments and concluded that the numerical damage could not be applied with better accuracy to the martensite phase because the indents are the same size or bigger than martensite grains. While there is a considerable difference in the indent and average grain size for ferrite, the identified parameters in their study could predict the damage in the ferrite phase [11].

Analytical models commonly applied to predict and plot FLDs by adopting major and minor limit strains are Swift, Chow–Hosford, Keeler–Brazier, and Marciniak– Kuczynski [11]. These approaches have been analyzed and compared in depth in research work by Basak and G Béres [12,13]. These approaches depend on some variables to develop the mathematical model and corresponding limit strains, i.e., Swift models require strain ratios and strain hardening coefficient values. Chow equations depend on the anisotropy coefficients, and for the Keeler–Brazier model, the thickness value and strain hardening coefficient are primarily important. The simulations in this study were performed using two approaches, i.e., the Keeler–Brazier, and M-K approaches, because of their efficiency, popularity, and simplicity in plotting FLDs of high-strength steels recommended by Duancheng Ma, Basak, and Béres. Kuang-Hua Chang et al. [14] used damage percentage detection to show the marginal zone of an FLD, where the area under this zone is safe, while the area above it is a damage zone.

Virtual modeling and simulation tools for the sheet metal processes play an essential role in predicting mechanical behavior and have become an essential and inevitable part of each industry [15,16]. These tools introduce full behavior prediction models of the metals, starting from their production to the heat treatment and testing processes. It helps engineers and designers to reduce the resource-consuming experiments and push the overall economic and technical aspects forward. A significant amount of work is being carried out toward developing, validating, and implementing numerical models to improve the accuracy of the simulation results. The numerical simulations method is an intelligent tool used to solve complex equations with different variables, called constitutive equations. These equations and variables represent the physical situations of the materials, e.g., deformation mechanisms, and help analyze the phenomena of mechanical deformation, damage, and failure. It can also predict crucial mechanical thresholds for high-end application materials during melting, casting, forming, and machining [17].

Hutchinson [18] built a model for FCC crystal, which was later expanded for BCC and HCP crystal structure and applied to DAMASK (Düsseldorf Advanced Materials Simulation Kit). Michel et al. [19] stated that the heterogeneity of ferrite and martensite phases in DP steels and their elastic stiffness coefficients play an essential role in the convergence behavior and stability schemes. After that, Diehl et al. [20] used some FFT assumptions and developed the stiffness coefficients applied later in the DAMASK framework. [7] The phenomenological crystal plasticity model is employed as plasticity law, which assumes that plastic deformation occurs on a slip system when the resolved shear stress exceeds a critical value. The resolved critical shear stress depends on the amount of the applied stress and can be calculated from a relation known as Schmid's law in Equation (1) [19,21].

$$\dot{\boldsymbol{\gamma}}^{\alpha} = \dot{\boldsymbol{\gamma}}\_{0} \left| \frac{\boldsymbol{\pi}^{\alpha}}{\mathbf{S}^{\alpha}} \right|^{n} \text{sgn}(\boldsymbol{\pi}^{\alpha}) \tag{1}$$

Analytical techniques can merely solve partial differential equations used in these constitutive models more straightforwardly, while numerical methods are essential for complex forms. Many efforts have been made to reach a framework connecting these boundary problems with physical phenomena [18]. Many numerical methods, i.e., Finite Element Method (FEM), Finite Volume Method (FVM), spectral method, and the Fast Fourier Transform (FFT)-based (Crystal Plasticity Finite Element Method) CPFEM method are usually used [22,23]. The difference between FE and SP methods lies in their homogenization technique and consequent time saving, as Shanthraj et al. [24] reported. The analytical FLD strain calculation models must be solved numerically with a high-efficiency solver to construct FLDs [25].

A representative volume element (RVE) was constructed for this purpose as a virtual sample to express the properties and microstructure of the material in the sample. In this study, spectral solver methods were applied using Fast Fourier Transform (FFT) to solve the boundary value problem for mechanical equilibrium and damage phase field. An FFT-based spectral solver shows higher efficiency and faster computational time over the domain.

The elastic and plastic material properties are governed by the crystal plasticity-based constitutive model, which is extensively used to study the deformation in crystalline materials. Besides considering microstructural parameters, it also considers some fitting parameters to compensate for the influence of some complex phenomena happening during the processing route of a specific grade of the material [25,26]. The simulations presented in this work were performed using DAMASK, which is available as free and open-source software [25]. It aims to simulate the material using crystal plasticity principles within a finite strain framework for continuum mechanical considerations by modeling the material point (Fourier point) inside the constructed mesh. The plasticity laws implemented in DAMASK are represented by isotropic plasticity, phenomenological crystal plasticity, or dislocation density-based crystal plasticity. The phenomenological crystal plasticity models were used in these simulations of polycrystal models. Furthermore, as deformation of polycrystalline aggregate strongly depends on the respective orientations of grains inside the lattice, consideration of Euler angles and rotation matrices was taken during the implementation of the model [20,25].

In this work, the establishment of a crystal plasticity-based approach to evaluate the forming limits of multi-phase steels by numerical simulation has been carried out successfully. A plastic instability approach and a high-efficiency numerical solver were employed to investigate multi-phase materials' formability by adopting major and minor strains in FLDs. Specifically, the microstructural attributes of martensite, i.e., grain size and phase fraction under different loading conditions, were studied. Furthermore, the damage initiation and propagation were determined by a CPFFT-based spectral solver using a phenomenological model for plastic deformation and Hooke's law for elasticity. The conclusions of this study can help the sheet metal industries as a guide for the designers to process material in an improved way.

#### **2. Methodology**

The numerical simulation modeling approach is a systematic process that uses welldefined pre-processing steps, running a set of simulations, post-processing, and visualization of the results, as shown in Figure 1. Firstly, RVEs were constructed by varying martensite volume fraction and martensite grain size by using DREAM.3D [27]. Next, the DAMASK adopts individual grains details to run numerical simulations using given load and geometry configuration files. Constitutive equations and microstructural attributes

for ferrite and martensite phases were adopted from Tasan's work [10]. Ductile damage and degradation parameters already implemented by Shanthraj et al. [24] were applied to the ferrite phase within a pre-developed model by Roters and Tasan et al. [20,25]. Finally, the values of local and global results were extracted after simulations and visualized by ParaView and other trend-plotting tools.

**Figure 1.** A simple flow chart showing the methodology adopted and data flow in this study.

#### *2.1. Data Collection*

The microstructural attribute values, e.g., martensite fractions and grain sizes, were taken from Tasan et al. [10], as shown in Table 1. The study was carried out for DP600 with the variation of martensite attributes (fractions 1.2%, martensite grain size from 1.0 to 4.3 μm) and a considerable ferrite grain size difference within 2.2 to 14.5 μm.

**Table 1.** The chemical composition and microstructural characteristics of DP 600 used in the current study. Reprinted with permission from Ref. [10] Copyright 2014 Elsevier.


#### *2.2. RVE Construction*

From DP steel values in Table 1, four RVE models (A–D) were constructed using DREAM.3D with variations in RVE dimensions, phase fractions, martensite grain size, and spatial distribution (refer to Figure 2). However, in all the five RVE models, the same hardness properties of martensite and ferrite grain size were considered, as shown in Table 2.

**Figure 2.** Top surfaces of 3D RVE models of the dual-phase steels, where the black particles represent the martensite phase inside the grey ferrite matrix.

**Table 2.** The projected RVEs used in the study with ferrite grain size 8.4 ± 6.1 μm for all the RVEs.


Suitable statistical distributions, cubic crystal structure, and ellipsoid grain shapes closely mimic the actual microstructure of DP Steel. In addition, different dimensions of synthetic volume were adopted to check the effect of RVE size on the global stress–strain curves and FLDs. Qayyum et al. [28] have shown a more detailed framework for the RVE generation using the DREAM.3D pipeline in their work. If interested, the readers are encouraged to refer to their work for further details.

#### *2.3. Pre-Processing Stage*

#### 2.3.1. Material Properties

The ferrite and martensite phases in DP steel have some common elastic–viscoplastic properties, which help build a more manageable material file framework, despite variations in their mechanical behavior and properties. The material parameters and damage values of both phases were adopted from Qayyum et al. [7], wherein already developed and validated models from the framework of DAMASK [25] were adopted. The elastic coefficients, initial and saturated shear resistances of slip systems and fitting parameters from the already published literature [10,29] were used as presented in Table 3. Regarding damage, the already developed models from Roters et al. [25] were incorporated and adopted in the material configuration files of respective RVEs. The critical plastic strain value crit for the ferrite phase was taken as 0.5 [10].

#### 2.3.2. Boundary/Loading Conditions

Different loading conditions were used in the simulation process with four strain states on each RVE model along y-directions while controlling x-directions and freeing the z-directions (refer to Figure 3).

The periodic boundary conditions were stated as:

In the uniaxial tension state:

$$\begin{array}{c} \dot{\mathbf{F}} = \begin{bmatrix} \* & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \* \end{bmatrix} \times 10^{-3} \begin{array}{c} \mathbf{s}^{-1} \\ \end{array} \qquad \qquad \begin{array}{c} \mathbf{\&} \\ \end{array} \qquad \qquad \begin{array}{c} \mathbf{P} = \begin{bmatrix} 0 & \* & \* \\ \* & \* & \* \\ \* & \* & 0 \end{bmatrix} \text{Pa} \tag{2}$$


**Figure 3.** Projected loading states to conclude the forming limit diagram.

In the plane strain state:

$$\dot{\mathbf{F}} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \* \end{bmatrix} \times 10^{-3} \text{ s}^{-1} \qquad\qquad\text{ @ }\qquad \mathbf{P} = \begin{bmatrix} \* & \* & \* \\ \* & \* & \* \\ \* & \* & 0 \end{bmatrix} \text{ Pa} \tag{3}$$

In the partial biaxial stretching state:

$$\begin{array}{c} \mathbf{F} = \begin{bmatrix} 0.5 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \* \end{bmatrix} \times 10^{-3} \begin{array}{c} \mathbf{s}^{-1} \\ \end{array} \qquad \qquad \begin{array}{c} \mathbf{P} = \begin{bmatrix} \* & \* & \* \\ \* & \* & \* \\ \* & \* & 0 \end{bmatrix} \text{Pa} \end{array} \tag{4}$$

In the equi-biaxial stretching state:

$$\begin{array}{c} \mathbf{F} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \* \end{bmatrix} \times 10^{-3} \begin{array}{c} \mathbf{s}^{-1} \\\\ \end{array} \qquad \qquad \begin{array}{c} \mathbf{P} = \begin{bmatrix} \* & \* & \* \\ \* & \* & \* \\ \* & \* & 0 \end{bmatrix} \text{Pa} \end{array} \tag{5}$$

where **F˙** and **P** are the rate of the deformation gradient and first Piola–Kirchhoff stress tensors, respectively. The coefficients denoted by '∗' express the stated complimentary conditions. Using these conditions, the four strain states were applied in the y-direction, with a 1 × 10−<sup>3</sup> s−<sup>1</sup> iso-static strain rate, as shown in Figure 3. For readers not familiar with the modeling strategy, a brief model background is provided in the work of Qayyum et al. [7] and Duancheng Ma et al. [30].

#### *2.4. Evaluation of FLDs by M-K and K-B Approaches*

#### 2.4.1. M-K Approach

In this approach, the engineering stress and strain values were extracted by applying customized subroutines. Then engineering stress–strain curves of the four loading conditions for each RVE were plotted. According to Duancheng Ma [30], Drucker's stability criterion claims that the forming limit occurs at the maximum stress (localized necking point) on the engineering stress–strain curve.

#### 2.4.2. Keeler-Brazier Approach

In this approach, the true stress and strain values were extracted from applying already customized subroutines, and then the true stress–strain curve was plotted. As per the Keeler–Brazier model detailed in appendix A, the major strain values (E1) depend on t (RVE thickness), n (strain hardening coefficient of material), and E<sup>2</sup> (minor strain values detected from the tensor matrix at an increment, corresponding to the maximum stress on the stress and strain curve). Therefore, the major and minor strain values were available for the four loading conditions. These and other extended details about post-processing and an activity flow chart of M-K and K-B approaches are given in Appendix A.

#### **3. Results**

Simulations of multi-phase DP steel were processed for 3D RVEs with damage consideration, wherein four models were virtually constructed to detect the variation in FLDs' behavior by varying microstructural attributes. By simulating different synthetic volumes of RVEs with the same ferrite grain sizes, the varying behavior of engineering stress–strain curves and FLDs was observed.

#### *3.1. Effect of Martensite Fractions on Stress–Strain Curve and FLDs*

The numerical simulations were carried out for models A and B as per Figure 2, and the difference in martensite volume fractions, as shown in Table 2. The behavior of stress–strain curves of uniaxial, plane strain, biaxial, and equi-biaxial loading is shown in Figure 4. By increasing martensite volume fractions, the yield and ultimate strengths increased while the values of corresponding engineering strains decreased. The stress values of uniaxial loading were almost like those of plane strain loading conditions but much lower than biaxial and equi-biaxial loading conditions. By calculating areas under engineering flow curves for the four loading conditions, the plastic work per unit volume increased by decreasing the martensite fraction in DP steels. The plastic work performed by the biaxial loading case ranged between 37 and 39 MJ/m3, while uniaxial and plane strain loading cases ranged from 28 to 36 MJ/m3.

**Figure 4.** Comparison of stress–strain curves of different martensite volume fractions in (**a**) uniaxial tension, (**b**) plane strain tension, (**c**) partial biaxial tension, and (**d**) equi-biaxial loading.

Concerning FLDs, the results are more significant, as shown in Figure 5. It was observed that the increase of martensite volume fractions decreased the values of major and minor strains of FLDs slightly in the M-K approach. However, in the case of the K-B approach, the formability difference was not observed.

**Figure 5.** The FLDs in martensite fraction differences in the case of (**a**) M-K approach and (**b**) K-B approach.

By plotting the logarithmic true stresses and strains of models (A and B), both variables were directly proportional to each other; when the values of n-values were detected from the slope of the graph. It was found that by decreasing the martensite volume fractions, the values of n-values increased, as shown in Figure 6 and Table 4. The trend is similar to what is already reported in the literature [25].

**Figure 6.** Comparison of strain hardening exponent (n) with corresponding martensite volume fraction, i.e., 18% and 17%.


#### *3.2. Effect of Martensite Grain Size on Stress–Strain Curve and FLDs*

By simulations of models (A, C, and D) as per Figure 2, the influence of the difference in martensite grain sizes was plotted as stress–strain curves of the four loading conditions in Figure 7. The stress–strain curves do not significantly affect the variation of martensite grain size in uniaxial loading. Contrarily, it was observed that the maximum stress of medium and small martensite grain size was almost the same, but the RVEs with bigger martensite sizes failed at lower strain values.

**Figure 7.** *Cont*.

**Figure 7.** Comparison of stress–strain curves of martensite grain sizes in (**a**) uniaxial tension, (**b**) plane strain tension, (**c**) partial biaxial tension, and (**d**) equi-biaxial loading.

By increasing martensite grain sizes, the limiting major and minor strains decreased in the case of the M-K approach with some exceptions in uniaxial loading cases, as shown in FLDs in Figure 8. On the other hand, in the case of the K-B approach, the variation in martensite grain sizes did not have a remarkable effect on formability.

**Figure 8.** The FLDs in martensite grain sizes in the case of (**a**) M-K approach and (**b**) K-B approach.

#### *3.3. Necking Band*

In the case of the M-K approach, as shown in Figure 9, the FLD band was observed for models A and B, where the green line is the necking start at 0% of material degradation, the blue line is localized necking point, and the red line represents the fracture of the RVE at 20% of material degradation. The area under the green line is a safe zone, and above the blue line is the damage zone, where the necking starts.

While in the case of K-B, by using models A and D, the values of major strains at the plane strain loadings were almost the same when the damage values were changed, while the values of the major and minor strains in other loading conditions changed. In addition, the marginal zone of the FLD was not observed for either model, as shown in Figure 10.

**Figure 9.** FLDs are showing marginal zones of (**a**) model A and (**b**) of model B.

**Figure 10.** The FLDs marginal zones for (**a**) martensite grain sizes = 2.7 μm and for (**b**) martensite grain sizes = 6 μm.

#### *3.4. Local Damage Evolution*

Damage initiation in models A and C were analyzed on the top surfaces of 3D RVEs, as shown in Figures 11 and 12, respectively. The frames were extracted during post-processing of the simulation results to compare local behaviors of the RVEs at the necking point, which showed 20% and 30% of global material degradation, respectively. The local damage field is presented with the help of values ranging from 1.0 to 0.0, i.e., undamaged to fully damaged state, respectively.

**Figure 11.** Damage behavior of dual-phase steels in model A with martensite grain sizes = 2.7 μm.

In model A (refer to Figure 11), the voids coalesced at 45 degrees of the load direction for uniaxial tension. The voids propagated sharply in this direction by the application of further load. In relation to plane strain loading, the damage initiation occurred at 45 and 90 degrees of load direction. With additional loading, the voids propagated as sharp straight lines at 90 degrees because of fixed loading at x-directions.

At partial biaxial loading, the damage was initiated by large voids. These voids propagated in 45 and 90 degrees to load directions. In contrast, in the equiaxial loading case, the damage propagated in different and random directions, such as 0, 45, and 90 degrees to load directions (refer to Figure 11). In the case of model C (refer to Figure 12), the local damage evolution did not differ much from model A, especially in the uniaxial case. However, it was observed in the plane strain case that the larger martensite grains constraint the damage behavior, and the damage tended to form around the martensite grains. The damage propagation also behaved in the same manner in the case of partial-biaxial and equi-biaxial loadings. The damage was initiated at the inter-granular level but continued growing trans-granularly with increased matrix degradation. In addition, stress relaxation near the damage areas affected the damage propagation and strain localization.

**Figure 12.** Damage behavior of dual-phase steels in model C with martensite grain sizes = 4 μm.
