1 Introduction

Remarkable need of the automotive and aerospace industries for new metallic materials that can meet stringent requirements regarding weight/property ratio has recently been observed. This need is a driving force for fast development of modern innovative steel grades. The most widely used example of those advances steels are dual-phase (DP) steels with the tensile strength of 400 to 1200 MPa. The name dual phase was first reported by Hayami and Furukawa in Reference 1. The DP steels have been successfully applied in the production of the automobile structural parts because they are characterized by a combination of good formability, high bake hardenability, and crash worthiness. These elevated properties are the results of sophisticated microstructure morphologies, which consists mainly of a ferrite matrix (around 70 to 90 pct) and a hard martensitic phase (around 10 to 30 pct). However, it has to be mentioned that small amounts of bainite, perlite, or retained austenite may also be present in the DP microstructure. Properties of DP steels are affected by many factors, including: volume fraction of martensite, average carbon content and carbon distribution in martensite, ductility of martensite, distribution of martensite, ferrite grain size, alloying elements content in ferrite etc.[2,3]

The most usual approaches to obtain a DP microstructure are controlled cooling after hot rolling (Figure 1(a)) or continuous annealing of cold-rolled sheets (Figure 1(b)). The former approach is usually used to produce thick DP steel strips, while the latter is used for manufacturing thin sheets often applied in the automotive industry.

Fig. 1
figure 1

(a) Hot rolling and controlled cooling, (b) cold rolling and continuous annealing[4]

Full size image

Nowadays, while designing new manufacturing technologies, engineers are considering not only global homogenous material behavior but also they try to incorporate micro scale phenomena like microstructure evolution or failure. As mentioned, due to strict market demands for final products with reduced weight and increased strength properties, the problem of material failure during manufacturing stages is becoming crucial. Zones where fracture can initiate during production stages should be identified and manufacturing cycle should be redesigned to avoid such behavior before industrial trials. Also extended service life of modern products require determination of failure probability during exploitation conditions. Experimental analysis can provide all the required information; however, it is time consuming and expensive at the same time. Thus, to reduce costs of development of manufacturing technology for new products, advantages of computer aided design are more often used in industrial practice.

2 State of the Art in the Numerical Modeling of Fracture in DP Steels

Modeling of fracture in DP steels is a complex task because of the composite character of the investigated microstructure that consists of two phases with significantly different mechanical properties. Multiphase microstructure in these steel grades results in difficulties during defining universal stress intensity factor (K IC) or J-integral factor. There are many experimental investigations where the identification problem of appropriate fracture toughness parameter for different DP steel grades was undertaken.[5–11] Numerical modeling with applied fracture toughness parameter based on the real experimental procedure gives usually satisfactory results. Such parameter can be used with numerical investigation only after extensive experimental investigation. Another approach is used when large deformation is induced into the material during production. This requires accurate definition of the forming limit diagram (FLD),[12] which gives information about critical deformation leading to fracture. Such criterion can be simply adopted into the numerical application but the FLD has to be determined separately for each DP steel grade under investigation. Examples of modeling fracture with the FLD approach are presented in (Reference 13 through 15). Other works[16,17] based on the maximum shear strain criterion take into account fracture generated during mechanical shearing process. Lack of accuracy of these models and discrepancies observed between model prediction and experimental measurements forced scientists to develop multi scale models with accurate description of microstructure morphology. Thus, the multi scale models based on the virtual material representation of microstructure are being intensively developed.[18–21]

There are three approaches to the virtual material description called representative volume elements (RVE), unit cell (UC), and digital material representation (DMR). The RVE is a model of the material that is used to determine the corresponding effective properties for the homogenized macroscopic model. The RVE should be large enough to contain sufficient information about the microstructure in order to be representative; however, it should be much smaller than the macroscopic body.[22 23]

The UC is a part of RVE that enables obtaining results for particular part of the material. Thus, the UC is not representative for the whole numerical model. Application of the UC idea, enables analyzing material behavior in particular location, e.g., crack initiation along the inclusion, without focusing on the rest of the material. As a result the UC provides data only for local analysis, not for the whole material (unlike the RVE). Usually several UCs can be considered as the RVE.

However, both presented approaches can have a detailed or simplified geometry of microstructure features. In the simplified model e.g., only similar volume fraction of particular phase is considered while the shape of this phase is not regarded. The model can provide representative global response, while morphology is significantly simplified and only statistically similar.[24–26] When detailed representation of morphology of real microstructure is needed, the term of the digital material representation is used. The DMR has already been introduced in many research works.[20,27–29] The definition according to[30] states that Digital Material is a material description based on measurable quantities that provides the necessary link between simulation and experiment.

So, the main objective of the DMR is to create digital representation of microstructure with its features represented explicitly to match real microstructure morphology. In that sense, the DMR can be used as UC or, if it meets specific criteria, it can be considered as the RVE (Figure 2).

Fig. 2
figure 2

Concept of the RVE, DMR, UC, and SSRVE

Full size image

The DMR concept was recently proposed and is dynamically evolving.[29–32] As mentioned, the main objective of the DMR is creation of the digital representation of microstructure with its features (i.e., grains, grain orientations, inclusions, cracks, different phases, etc.) represented explicitly. Generation of material microstructure with specific geometrical features and properties is one of the most important algorithmic parts of the DMR approaches. The DMR is further used in numerical simulations of processing or simulation of behavior under exploitation conditions and the more accurate, in the case of geometry and properties, the digital representation is the more accurate results can be obtained. Due to that, a lot of research is directed towards development of methods responsible for creation of 2D and 3D representations of analyzed microstructures.

To obtain an accurate description of the 2D microstructure an image processing (IP) methods are usually applied. As input data for this analysis an SEM/EBSD image can be used. That way not only information regarding microstructure geometry is obtained but also information about initial crystallographic orientation is provided. This approach was successfully used in References 33, 34. Unfortunately, the approach is time consuming and expensive, because each numerical simulation based on the DMR require SEM/EBSD analysis. That is why IP is also applied to the optical microscopy images that are much more affordable. However, in this case, only information regarding grain morphology is obtained see e.g.[35] Presented approaches are even more demanding when 3D digital representations are required. In 3D cases, the DMR is usually created based on the reconstructed 2D slices obtained using a destructive method. Again an optical or scanning electron microscopy can be used during the serial sectioning procedure to provide input data for IP and reconstruction algorithms. The conventional approach to serial sectioning based on the manual labor is extremely time consuming and requires a series of subsequent polishing and optical or scanning electron microscopy operations.[36] The solution may be application of specially designed equipment e.g., Robo-Met.3DTM[37,38] that automatically provides stack of 2D images that are subjected to 3D reconstruction algorithms. This equipment is limited only to optical microscopy images. To obtain a series of 2D images from the scanning microscopy in an automated manner a SEM/FIB/EBSD technique can be used. The focused-ion beam (FIB) plays a crucial role in this procedure. The major advantage of the method is possibility to obtain not only 3D morphology of particular grains but also information about their crystallographic orientation.[39] Unfortunately, relatively small areas of material can be reconstructed by this approach.

Methods presented above provide an exact 2D or 3D digital representation of analyzed microstructures; however, they require a series of experimental research and metallographic analysis. Thus, a series of numerical methods based on the cellular automata, Monte Carlo or Voronoi Tesselation was also developed.[31,32,40]

The IP methods were used in the present work to develop the DMR model that can be used to explicitly predict influence of the martensite fraction on fracture behavior.

3 Digital Material Representation Based on Optical Microscopy Images

In order to exactly recreate microstructure morphology of the DP steel (Table I), a series of heat treatment operations was conducted to provide samples with different martensite fraction for further metallographic analysis and IP.

Table I Chemical Composition of the Investigated Steel
Full size table

Four various quenching conditions were applied to differentiate the amount of martensite namely: 1, 52, 105, and 475 °C/s. As a result, four different microstructure morphologies were obtained with martensite volume fractions: 27.7, 43.3, 57.7, and 64.5 pct, respectively, as seen in Figure 3. The morphology of martensite is closely related to e.g., cooling conditions, state of the initial microstructure, level of alloying elements, etc. Thus, the problem of representativeness of investigated microstructures both from experimental and numerical point of view remains open.[41] That is why for the purpose of development of the DMR models Authors selected microstructures located in the center of the plate where the most uniform processing conditions are expected.

Fig. 3
figure 3

DP steel microstructures after various cooling rates

Full size image

The developed IP algorithm was used to transfer obtained optical microscopy images into the DMR form. First, microstructure image (Figure 4(a)) is subjected to digital treatment with the threshold function. This step is realized in the GIMP software (GNU Image Manipulation).[42] The thresholding algorithm performs simple binarization procedure. White pixels represent the pixels of the image, which values are below the threshold range, and black pixels represent pixels with values higher than the threshold range. After thresholding, some noise can be observed in the image, which disrupts visual separation of two phases (Figure 4(b)). To remove the noise, the filtering algorithm has to be applied (Figure 4(c)). Finally erosion algorithm is used to remove grain boundaries from the microstructure and leave only martensite islands in ferritic matrix as seen in Figure 4(d).

Fig. 4
figure 4

(a) Real microstructure image, and DMR images after (b) thresholding (c) filtering, and (d) removing grain boundaries algorithms

Full size image

The binary form of the microstructure presented in Figure 4(c), is an input for the second stage of the algorithm—separation and coloring procedure to identify subsequent ferrite grains. For this purpose, the cellular automata (CA)[43] based algorithm was used within the work. The algorithm involves coloring stage combined with the grain growth model. A single pixel in each grain is selected to represent grain nuclei. Next, a simple transition rule is applied: when a neighbor of a particular cell in the previous time step is in the state ‘already grown’, then this particular cell can also change its state. The grains grow with no restrictions until they fill the entire investigated domain (Figure 5).

Fig. 5
figure 5

Image coloring procedure

Full size image

As a result, each grain has the unique color identifier, what is shown in Figure 6.

Fig. 6
figure 6

Image after coloring procedure

Full size image

Another stage is required to remove thin grain boundaries and is also based on the grain growth CA model. The developed method removes all black cells from the image. The transition rule states that when one or more of the neighbors surrounding the cell is black than the cell accepts a color of one of its neighbors (Figure 7).

Fig. 7
figure 7

Illustration of the grain growth CA model

Full size image
Fig. 8
figure 8

Effect of the CA grain growth algorithm

Full size image

The procedure is repeated until all black cells change their color. It has to be emphasized that in the present work different colors of subsequent grains do not represent crystallographic orientations, they are just used to distinguish particular grains. The outcome of the model is presented in Figure 8.

As a result, a temporary DMR without any black regions representing martensite islands is created. The last step of the algorithm is responsible for combining results obtained after simple IP (Figure 4(d)) with the temporary DMR from Figure 8. Eventually, complex DP microstructure morphology is restored, where both martensite islands as well as subsequent ferrite grains are clearly visible (Figure 9).

Fig. 9
figure 9

Main steps of the developed image processing approach

Full size image

Due to the fact, that only part of the entire macroscopic sample is taken into account in the simulations, in order to ensure the continuity of the design space, the periodic boundary conditions (PBC) have to be introduced. Unfortunately, image-based digital microstructures do not have periodic character. Thus, a buffer zone approach[44] that permits to induce periodicity has to be used as presented in Figure 10.

Fig. 10
figure 10

The DMR with the buffer zone

Full size image

The DMR is then an input for the mesh generation algorithm that performs discretization of the calculation domain. For this stage, the in-house code for finite element non-uniform mesh generation was developed. Details of the DMRmesh software can be found in References 34, 45, 46, and only major steps are presented below:

  1. (a)

    Establishing border nodes on the basis of the DMR morphology. These nodes map the boundary between phases or grains.

  2. (b)

    Generation of nodes in the interior of investigated microstructure features. Each node is described by a defined radius and within this range no other node can be added. Additionally, nodes can have different radii and that way, a specific node distribution along the feature boundary can be obtained.

  3. (c)

    Triangulation with the Delaunay algorithm to create non-uniform meshes.

  4. (d)

    Improvement of the mesh quality by the Laplace smoothing and edge-swapping approaches.

Finally conforming mesh dedicated to numerical simulation that exactly replicates grain geometry of DP microstructure can be generated as seen in Figure 11. The DMRmesh software allow to generate triangular as well as quadratic finite element meshes.

Fig. 11
figure 11

Finite element mesh refined along phase/grain boundaries generated on the basis of DMR

Full size image

In case of PBC, additionally there must be the same number of nodes on both sides of the sample. Nodes from one side of the material are then connected with nodes on the opposite side. Thus, there is a two stage procedure:

  • First, the FE mesh is generated that ensures matching nodes positions on corresponding edges of the sample (Figure 12).

    Fig. 12
    figure 12

    Illustration of the PBC applied to the RVE: the west-east version

    Full size image

  • Second, a set of constrains is automatically applied to each pair of nodes by a developed Python script. The Python script was used due to a large number of constrains needed to be imposed. Each node along the corresponding edges is described by equations calculating displacement values at vertical edges U = [U x, U y]:

    $$ U_{x,y}^{\text{West}} = U_{x,y}^{\text{East}} + U_{x,y}^{\text{SE}} , $$
    (1)

    where U West,East x,y —displacements calculated in the particular nodes, U SE x,y —transfers displacement values between corresponding nodes.

As seen in Figure 12, only west-east PBC were applied in the paper to better reflect tensile conditions.

For all the simulations, the same material properties for the martensite and ferrite phases were used. At this stage of the research, the same material properties were assigned to the phases in subsequent DMR models. Such simplification can be done because main target was focused on the investigation the influence of the martensite volume fraction on failure. Based on simple rule of mixture, different material definitions were adopted to the surrounding buffer zone for the investigated cases. As a result, material properties of the buffer zone are associated with percentage amount of the ferrite and martensite and can be considered as properties of homogenous material. Number of finite elements used for discretization was set to approx. 500,000 to minimize the mesh sensitivity effect. Four node bilinear plane strain quadrilateral finite elements (CPE4R) were chosen for the discretization purposes. The applied specific mesh was obtained with mentioned FE mesh generation software.[25] Created DMR models are presented in Figure 13.

Fig. 13
figure 13

DMR models of DP steels obtained after cooling with rates: (a) 1 °C/s, (b) 52 °C/s, (c) 105 °C/s, and (d) 475 °C/s

Full size image

DMR models were created in commercial ABAQUS application and calculated using explicit solver. Boundary conditions presented schematically in Figure 14 were applied during calculations.

Fig. 14
figure 14

Boundary conditions used during calculations

Full size image

A combination of the above described steps leads to creation of DMR models incorporated into the commercial finite element software for further modeling of fracture of DP steels during deformation. An important step in creation of described numerical models, is assignment of the material properties to the ferrite, martensite as well as to buffer zone, respectively. The material model for the present investigation was identified on the basis of micropillar compression tests.

4 Identification of DMR Material Models

In the last decade, micropillar compression test has been developed to characterize mechanical properties of materials at the micro/submicron scale.[47–49] It is a miniaturized version of conventional uniaxial compression test that can determine mechanical behavior of nano-/micro-structured materials by compressing micropillars under uniaxial condition. In the present research, micropillar compression tests were carried out to identify mechanical behavior of ferrite and martensite grains separately. Micropillars of ~1 μm in diameter were fabricated by using FIB in a dual beam scanning electron microscope (SEM, FEI Quanta 200). A top-down methodology was adopted to mill micropillars using a Ga+ ion beam operated at 30 kV. The beam current was reduced from 5 nA for coarse milling down to 50 pA for final polishing, which was designed to minimize the surface damage caused by the Ga+ ion beam. Then micro compression tests were carried out by using a diamond flat punch (5 μm in diameter) within a nanoindentation apparatus (UMIS, CSIRO). The flat punch was driven in a force-controlled mode to compress the pillars, and the displacement rate was ~2.0 nm/s. Figure 15 shows micropillars of investigated DP 600 steel cut out from single ferrite and martensite grains before and after deformation. The DP 600 is a very popular DP steel grade used in automotive industry mainly for front longitudinal elements in the car that can be subjected to severe deformation.

Fig. 15
figure 15

FIB-prepared micropillars of DP steels. FP means ferrite pillars, and MP means martensite pillars

Full size image

The loading force P and total displacement u total can be directly obtained from the micro compression tests. The initial length and cross-sectional area of the pillars can be measured from SEM images. As a result, the stress–strain relationship necessary for the finite element calculations can be obtained. The measured and approximated true stress–strain curves recorded during pillars deformation are shown in Figure 16.

Fig. 16
figure 16

True stress-true strain curves for DP steel pillars

Full size image

Obtained material properties were assigned to subsequent DMR models investigated within the paper.

5 Numerical Results

Eight numerical models were created on the basis of four investigated DMR models of DP steel with different levels of martensite volume fractions. The first set of four models did not take into account influence of fracture. The second set, incorporated Johnson–Cook ductile crack model described in (Reference 50, 51). In the ductile criterion, it is assumed that the equivalent plastic strain at the onset of damage is a function of stress triaxiality and strain rate:

$$ \varepsilon_{\text{iD}}^{\text{pl}} \left( {\eta ,\dot{\varepsilon }_{\text{i}}^{\text{pl}} } \right) = \frac{{\varepsilon_{\text{i}}^{ + } { \sin }h\left[ {k_{0} \left( {\eta^{ - } + \eta } \right)} \right] + \varepsilon_{\text{i}}^{ - } { \sinh }[k_{0} (\eta - \eta^{ + } )]}}{{{ \sin }h[k_{0} (\eta^{ - } + \eta^{ + } )]}}, $$
(2)

where ɛ + i and ɛ −i —equivalent fracture strains for equibiaxial tensile and equibiaxial compressive deformation, respectively, η—stress triaxiality (a ratio of the equivalent mean stress σ m to the Misses equivalent stress σ i), k 0—parameter obtained experimentally. These parameters depend on the material, strain rate, and temperature of the process. Failure in this model occurs when state variable defined by Eq. [3] reaches value of 1.

$$ w_{\text{D}} = \mathop \smallint \nolimits \frac{{{\text{d}}\varepsilon_{\text{iD}}^{\text{pl}} }}{{\varepsilon_{\text{iD}}^{\text{pl}} \left( {\eta ,\dot{\varepsilon }_{\text{i}}^{\text{pl}} } \right)}} = 1, $$
(3)

where \( {\text{d}}\varepsilon_{\text{iD}}^{\text{pl}} \)—plastic strain increment per simulation unit time.

The crack initiation parameter for the ferrite failure was adopted in the present work based on[51] and its evolution is presented in Figure 17.

Fig. 17
figure 17

Derived damage curve for the ferrite in DP steel[51]

Full size image

Examples of obtained results based on the DMR approach without failure model and combined with ductile failure criterion are presented in Figures 18 and 19, respectively.

Fig. 18
figure 18

Strain distributions without fracture model for microstructure obtained after: (a) 1 °C/s, (b) 52 °C/s, (c) 105 °C/s, and (d) 475 °C/s cooling rates

Full size image
Fig. 19
figure 19

Strain distributions with the fracture model for microstructure obtained after: (a) 1 °C/s, (b) 52 °C/s, (c) 105 °C/s, and (d) 475 °C/s cooling rates

Full size image

Because real microstructure morphology is taken into account highly inhomogeneous strain distribution is obtained during deformation. Major part of accumulated plastic strain is within the softer ferrite phase. Zones where strain value is high e.g., due to morphological features can be the locations of ductile crack initiation. In models with different martensite volume fractions, level of the strain localization during deformation between martensite islands is significantly dissimilar. As a result, in the model with lower martensitic volume fraction cracks start to initiate and propagate in the end of the deformation. On the other hand, in the model with high martensite volume fraction due to higher strain localization resulting from large martensite area fraction, cracks start to initiate earlier during the deformation. Presented comparison clearly highlights how it is important to consider failure models in simulations of such complex microstructures.

However, not only detailed information regarding inhomogeneous strain distribution, and fracture propagation during deformation can be obtained. Using this approach, a response in the form of flow stress model by homogenization can be obtained as seen in Figure 20. The main idea of the homogenization is to relate the micro fields to the macroscopic constitutive equations:

$$ \begin{gathered} S\left( x \right) = \frac{1}{V}\mathop \smallint \limits_{v}^{ } \sigma \left( x \right){\text{d}}V \hfill \\ E\left( x \right) = \frac{1}{V}\mathop \smallint \limits_{\text{v}}^{ } \varepsilon \left( x \right){\text{d}}V ,\hfill \\ \end{gathered} $$
(4)

where σ—stress value in each element, ɛ—strain value in each element, V—the total area occupied by a given phase.

Fig. 20
figure 20

Average stress–strain charts with and without fracture included

Full size image

Results presented in Figure 20 confirm in a quantitative manner results from Figures 18 and 19. Future work should also take into account the brittle failure mode in martensite, which can also influence presented behavior. It has to be emphasized that the major advantage of data presented in Figure 20 lies in the fact, that these flow curve models can also be used during numerical simulation of DP steel behavior at the macro scale level.

6 Conclusions and Future Work

Inhomogeneous strain distribution and fracture propagation in the microstructure of DPs steel simulated on the basis of DMR were presented in the paper. The major conclusion from the research is that it is important to include representation of geometry of the phases during numerical simulations when complex AHSS are investigated. The IP combined with the CA grain growth model provides a possibility to accurately transfer morphology into the digital form. Precise consideration of shape of phases received after quenching process allow to predict inhomogeneous strain distribution across microstructure during further processing operation. The strain localization zones seems to be the locations for ductile cracks initiation. It can also be concluded that by application of the DMR with the ductile criterion, it is possible not only to simulate crack initiation and subsequent propagation but also make relations between micro scale behavior and macroscopic response. Presented results clearly show that the martensite has noticeable influence on the material fracture behavior and it must be taken into account during numerical modeling of DP steel grades.