Studying the Damage Evolution and the Micro-Mechanical Response of X8CrMnNi16-6-6 TRIP Steel Matrix and 10% Zirconia Particle Composite Using a Calibrated Physics and Crystal-Plasticity-Based Numerical Simulation Model

Abstract

The mechanical behavior of newly developed composite materials is dependent on several underlying microstructural phenomena. In this research, a periodic 2D geometry of cast X8CrMnNi16-6-6 steel and 10% zirconia composite is virtually constructed by adopting microstructural attributes from literature. A physics-based crystal plasticity model with ductile damage criterion is used for defining the austenitic steel matrix. The zirconia particles are assigned elastic material model with brittle damage criterion. Monotonic quasi-static tensile load is applied up to 17% of total strain. The simulation results are analyzed to extract the global and local deformation, transformation, and damage behavior of the material. The comprehensively constructed simulation model yields the interdependence of the underlaying microstructural deformation phenomena. The local results are further analyzed based on the interlocked and free regions to establish the influence of zirconia particles on micro-mechanical deformation and damage in the metastable austenite matrix. The trends and patterns of local strain and damage predicted by the simulation model results match the previously carried out in-situ tensile tests on similar materials.

1. Introduction

The need for developing materials with good strength to weight ratios which possess superior formability and mechanical properties is rising. These materials are in demand for automobile, construction, aviation, and several engineering applications [1,2]. Particle-based metal matrix composites have evolved as an alternative to complex alloys of conventional materials and provide various possibilities for several targeted applications [3].

Transformation induced plasticity (TRIP) steel matrix reinforced composites with magnesium partially stabilized zirconia (Mg-PSZ) particles are a relatively newer class of materials recently developed by the Institute of Materials Science, TU Bergakademie Freiberg [4,5]. They demonstrate superior work hardening and ductility during deformation due to the interplay of several phenomena [6,7]. The metastable $\gamma$-austenite due to applied strain transforms to $\acute{\alpha }$-martensite phase by passing through a highly deformed $\epsilon$-martensite phase [8,9]. This transformation results in local plastic deformation but later restricts the flow of dislocations by acting as hard second phase particles. The magnesium partially stabilized zirconia (Mg-PSZ) particles undergo transformation from a monoclinic to cubic phase due to applied external load and hence undergo slight plastic deformation [10,11]. Upon further increasing of applied external load, ductile damage initiates at the ceramic/matrix interface and around the $\acute{\alpha }$-martensite particles [12]. The in-situ test results of this material have revealed that the brittle damage of the zirconia particles also starts to occur, facilitating the further degradation and eventual rupture of the material [13,14].

Dislocation glide, twinning, transformation, and grain boundary motion (depending on the grain size and geometry) facilitate the material flow. Conversely, the transformation of ductile $\gamma$-austenite to hard $\acute{\alpha }$-martensite, dislocation pinning around the martensitic islands and on the ceramic/matrix interface restrict the material flow. This is caused by a reduction in the mean free path and results in work hardening of the material [9]. It is essential to comprehensively understand the interplay of these microstructural attributes so that they can be adjusted to obtain desired mechanical properties.

Much experimental work has already been carried out in the recent past to understand the fundamental underlying deformation mechanisms in TRIP Steel MMCs [15,16]. It is reported that the local crystallographic orientation and physical attributes of each phase play an essential role in plastic deformation evolution and therefore are important to consider for understanding the material behavior precisely [17]. However, with so many microstructural parameters to adjust, it becomes quite cumbersome and expensive to optimize the material for desired applications by using conventional experimentally driven manufacturing, testing, and characterization routes.

The development of an accurate microstructurally informed numerical simulation model could potentially provide a good solution for inexpensive and fast material development. Some work has been carried out to empirically model the average mechanical behavior [15,16,18,19] and information was extrapolated. Due to the absence of microstructural information, these methodologies are intrinsically limited and depend heavily on the experimental test data. This leads scientists and researchers back to the same challenge which is needed to be addressed in the first place.

The recent developments in crystal plasticity-based modeling techniques have shown the potential to numerically model such complex material behaviors [20,21,22]. DAMASK is developed as a flexible crystal plasticity-based phase-field modeling toolbox for numerical simulations of such microstructurally complex composites [23,24,25]. Wong et al. [26] and Madivala et al. [27] developed a physics-based TRIP/TWIP model and incorporated it into the DAMASK framework. The model was adopted and calibrated for modeling the global and local deformation behavior of 16-6-6 and 16-7-6 TRIP/TWIP Steels [28,29]. When comparing the simulation results with experimental observations, it was shown by several researchers that the model can predict material behavior in low strain regimes, while it over-predicts local stress and strain evolution at high strains [30,31].

One reason for the over-estimation of results is the absence of the appropriate and calibrated damage criterion for a specific phase [32,33]. Recently Shanthraj et al. [34], following Miehe et al. [30], proposed a simple ductile damage model for the phase-field approach and incorporated it in DAMASK. In addition, a brittle damage model [35] following Kim et al. [36] was also proposed and incorporated in the DAMASK framework for initial level implementation. Several researchers adopted these damage models separately, i.e., for the ductile matrix [37,38,39,40] and brittle ceramic particles [24] to highlight their possible use by comparing with the experimental observations. The fitting parameters of these damage models are calibrated by comparing the global flow curve with tensile test results.

For TRIP Steel MMCs, it has been shown that excessive strain hardening takes place in the metastable austenitic matrix during deformation. These materials have been widely studied under low strain conditions, and the influence of several micromechanical behaviors has been precisely documented [4,17]. It is reported that with increasing strain, the rate of transformation and mean free path for dislocations to move drops exponentially. At high strains, three kinds of damage nucleation phenomena start to take place:

• Ductile damage in matrix
• Ceramic/matrix interface decohesion
• Brittle cracking of the ceramic particles

However, there is no order of initiation between these competing stiffness degradation phenomena, and they are reported to evolve simultaneously. With increased deformation, void coalescence occurs, which results in the final rupture of the component.

The localization of these damage nucleation phenomena and their evolution are vastly interconnected with microstructural architecture and evolving physical attributes but have not been studied in detail before. The limitations of the experimental setups hinder such a study and limit the scope of material microstructure engineering for the TRIP Steel Mg-PSZ MMCs. A robust phase-field crystal plasticity-based numerical simulation model incorporated with calibrated ductile and brittle damage criteria can be used to analyze the local deformation and damage behavior of such materials. Such a model will also help in understanding the interaction of several micromechanical phenomena at high strains.

In the current work, a capability to address this challenging problem of modeling deformation and damage in such materials is presented. In Section 2, the construction of the RVE and numerical simulation methodology with model parameters is presented. For simplicity, a 2D RVE is loaded under a quasi-static uniaxial tensile load, and results are recorded. The results are presented in Section 3, showing the global and local behavior of the material in detail. The local results are further dissected to compare the difference of attributes in the interlocked and free matrix. Discussion of obtained trends and results compared with previous publications is comprehensively carried out in Section 4 of the article. Finally, the study is concluded in Section 5 with a brief outlook.

2. Materials and Methods

A periodic virtual 2D RVE was generated by adopting microstructural attributes from the literature [16]. For the X8CrMnNi16-6-6 austenitic steel matrix, the calibrated physics-based TRIP/TWIP model from the already published work of Qayyum et al. [28] was adopted with incorporated ductile damage criterion. For Mg-PSZ particles, a phenomenological model with elastic attributes was adopted [41] with an integrated brittle damage criterion. The parameters of both damage criteria were adjusted based on the critical plastic strain and critical strain energy for austenite matrix and zirconia particles reported in the literature [13], respectively. The simulations were performed by applying a monotonic quasi-static tensile load. The data was extracted for each frame after numerical simulation by using available subroutines. The results were post processed to obtain detailed global and local evolution of physical parameters to understand their influence on deformation and damage.

2.1. RVE Generation

MMC composition is considered with 10% of Mg-PSZ particles for the present case. The material manufacturing by sintering and characterization details are published elsewhere [16]. In the current work, the microstructural attributes are adapted from previously published work for virtually constructing the representative volume element using already published methodology [42].

The microstructural statistical details are used as input in the Dream3D [43] pipeline, and virtual grain size distribution data for each phase is generated using the overall homogeneous orientation distribution function (ODF). Afterwards, 90% austenite phase and 10% zirconia phase with equiaxed grain size distribution are defined. A 3D periodic RVE is constructed and later sliced to get a 2D layer for the simulations from the input data. The current RVE is of 100 µm × 100 µm dimensions with 1 µm² element size, which means that it has 10,000 individual solution points.

It has been shown by Qayyum et al. [42] earlier that the 2D RVEs show higher contrast of local property evolution, and the quantitative results can be slightly misleading, especially at higher strains. 2D RVE is adopted in the present work to keep the computation costs low and simplify the results for easier understanding and better presentation. The model will be scaled up for 3D RVE simulations in the future. Some grains of austenite matrix and zirconia particles in the RVE seem smaller than 10 µm. It is because due to slicing, the third dimension is lost. For such full phase simulations, it is an inevitable assumption.

The IPF map shown in Figure 1a shows the overall neutral/unbiased orientation distribution of the RVE considered for the current case. The second phase distribution is heterogeneous, which closely resembles the observed material microstructure [19]. There are areas where the austenitic matrix is interlocked between the zirconia particle clustering, and there are some areas where the austenite matrix is observed to be free/open. These zones were intentionally kept this way to analyze and report the effect of heterogeneous microstructural attributes on the local deformation evolution (more about this in the local results section).

2.2. Material Modeling

2.2.1. TRIP/TWIP Model with Ductile Damage

A dislocation density-based crystal plasticity model incorporating both transformation-induced plasticity (TRIP) and twinning-induced plasticity (TWIP) was proposed by Wong and Madivala et al. [26,27], respectively. In this model, the total deformation gradient $\mathit{F}$ is decomposed in two elastic ${\mathit{F}}_e$ and plastic ${\mathit{F}}_p$ parts. The total velocity gradient is defined as:

$$\mathit{L}=\dot{\mathit{F}}{\mathit{F}}^{-\mathbf{1}}$$

(1)

Following Kalidindi’s approach [44] with pure dislocation slip, contributions from mechanical twinning and phase transformation are accounted for in the plastic velocity gradient determination:

$${\mathit{L}}_p=\left(1-{\displaystyle \sum }_{\beta =1}^{N_{tw}}{f}^{\beta }-{\displaystyle \sum }_{\chi =1}^{N_{tr}}{f}^{\chi }\right){\displaystyle \sum }_{\alpha =1}^{N_s}{\dot{\gamma }}^{\alpha }{\mathit{m}}^{\alpha }\otimes {\mathit{n}}^{\alpha }+{\displaystyle \sum }_{\beta =1}^{N_{tw}}{\dot{\gamma }}^{\beta }{\mathit{m}}_{tw}^{\beta }\otimes {\mathit{n}}_{tw}^{\beta }+{\displaystyle \sum }_{\chi =1}^{N_{\mathrm{tr}}}{\dot{\gamma }}^{\chi }{\mathit{m}}_{tr}^{\chi }\otimes {\mathit{n}}_{tr}^{\chi }$$

(2)

$f^{\beta }$ is the twin volume fraction and $f^{\chi }$ is the transformed martensite volume fraction. $\mathit{m}$ and $\mathit{n}$ vectors represent the directions/plane normals on which shear occurs at a rate of $\dot{\gamma }$. The present plastic velocity gradient ${\mathit{L}}_p$ does not account for the subsequent slip within the twined or transformed regions. The shear rates for each phenomenon are calculated separately in the model.

The microstructure is divided into edge dislocation density ${\rho }_e$, dipole dislocation density ${\rho }_d$, the twin volume fraction $f_{tw}$, and the martensite volume fraction $f_{tr}$. To calculate the evolution of dislocation densities, the Orowan model [45] gives shear rate on the slip system $\alpha$ as:

$${\dot{\gamma }}^{\alpha }={\rho }_e{b}_s{v}_0\exp \left[-\frac{Q_s}{k_{b}T}{\left\{1-{\left(\frac{\left|{\tau }_{eff}^{\alpha }\right|}{{\tau }_{sol}}\right)}^p\right\}}^q\right]\mathrm{sign}\left({\tau }^{\alpha }\right)$$

(3)

where ${\tau }_{eff}^{\alpha }$ is resolved shear stress on a slip system $\alpha$, ${\tau }_{sol}$ is solid solution strength, $b_s$ is the Burgers vector, $v_0$ is the dislocation glide velocity, $Q_s$ is the activation energy for a specific slip system, $k_b$ is Boltzmann constant, $p$ and $q$ are the fitting parameters. The evolution of edge dislocation density is given as:

$${\dot{\rho }}_e^{\alpha }=\frac{\left|{\dot{\gamma }}^{\alpha }\right|}{b_s{\Lambda }_s^{\alpha }}-\frac{2{\widehat{d}}^{\alpha }}{b_s}{\rho }_e^{\alpha }\left|{\dot{\gamma }}^{\alpha }\right|-\frac{2{\check{d}}^{\alpha }}{b_s}{\rho }_e^{\alpha }\left|{\dot{\gamma }}^{\alpha }\right|$$

(4)

The evolution of dipole dislocation density is given as:

$${\dot{\rho }}_d^{\alpha }=\frac{2{\widehat{d}}^{\alpha }}{b_s}{\rho }_e^{\alpha }\left|{\dot{\gamma }}^{\alpha }\right|-\frac{2{\check{d}}^{\alpha }}{b_s}{\rho }_d^{\alpha }\left|{\dot{\gamma }}^{\alpha }\right|-{\rho }_d^{\alpha }\frac{4v_{climb}}{\left({\widehat{d}}^{\alpha }-{\check{d}}^{\alpha }\right)}$$

(5)

where the dislocation climb velocity can be calculated using:

$$v_{climb}=\frac{3G{D}_0\mathsf{\Omega }}{2\pi k_{b}T}\frac{1}{\left({\widehat{d}}^{\alpha }-{\check{d}}^{\alpha }\right)}\exp \left(-\frac{Q_c}{k_{b}T}\right)$$

(6)

Here, $D_0$ is the self-diffusion coefficient for fcc Fe, $\mathsf{\Omega }$ is the activation volume for the climb, and $Q_c$ is the activation energy for the climb. The maximum value of glide plane distance two dislocations can have to form a dipole is ${\widehat{d}}^{\alpha }$ and minimum distance required for the two edge dislocations to annihilate ${\check{d}}^{\alpha }$ are calculated respectively as:

$${\widehat{d}}^{\alpha }=\frac{3G{b}_s}{16\pi \left|{\tau }^{\alpha }\right|}$$

(7)

$${\check{d}}^{\alpha }=C_{anni}{b}_s$$

(8)

where, $C_{anni}$ is a fitting parameter. The current model is based on the assumption that the twin formation occurs by the passage of $\frac{1}{6}11\overline{2}$ Shockley partials on every $\left\{111\right\}$ plane. Whereas the formation of martensite (fcc to hcp transformation) is obtained by the passage of $\frac{1}{6}11\overline{2}$ Shockley partials on every second $\left\{111\right\}$ plane. These are attributed as glide systems of higher-order as they require more energy for activation. Part of the model presented above has been taken from the previously published work of Wong and Madivala et al. [26,27]. The readers are encouraged to refer to their work for more detail about the model and assumptions.

In the work of Wong and Madivala et al. [26,27], the model was calibrated for high manganese TRIP/TWIP Steels, which possess high SFE. The phase transformation in such materials is γ → ε martensite transformation. However, in low SFE austenitic steels (currently under discussion), γ → ε → ά martensitic transformation occurs [8]. It is reported that there is a slight difference in the nano hardness of ε and ά martensite, as well as the crystallographic structure of these phases [46]. However, considering the problem’s complexity and model limitations, both ɛ-martensite and α´-martensite are treated as a single-phase for simplicity. With this assumption, the previously published TRIP/TWIP model [28] adjusted the physical parameters as per cast X8CrMnNi16-6-6 TRIP Steel, and the fitting parameters of the twinning, transformation, and dislocation glide are adopted. The physical and fitting parameters are shown in

Table 1. Optimized constitutive model parameters for X8CrMnNi16-6-6 TRIP steel. The table has been adapted from a previous publication [28] and is being reprinted with the permission of MDPI.

	Symbol	Description	Value	Unit
Martensite Transformation Parameters	Cthreshhold-trans.	Adj. parameter controlling trans threshold stress	0.5	--
	Ltr	Width of martensite lath during nucleation	5.0 × 10−8	m
	ttr	Average martensite thickness	5.0 × 10−6	m	[47]
	ß	ß-exponent in transformation formation probability	3.0	--
	CMFP-trans.	Adj. parameter controlling trans mean free path	10	--
Twinning FormationParameters		r-exponent in twin formation probability	1.0	--
	CMFP-TWIN.	Adj. parameter controlling twin mean free path	5.0	--
	Cthreshhold-TWIN.	Adj. parameter controlling twin threshold stress	1.3	--
		Stacking fault energy	10	mJ/m²	[48]
		Change in Gibbs free energy	−2.54 × 107	J/m3	[49]
Dislocation Glide Parameters		Solid solution strength	5 × 107	Pa
		Initial dislocation density	1.0 × 1012	m−2	[26]
		Adj. parameter controlling dislocation mean free path	55	--
		Avg. grain size	50 × 10−6	m	[50]

. The elasticity matrix constants for the austenitic phase and the transformed martensitic phase are provided in

Table 2. Single-crystal elastic constants for Austenite, Martensite, and Zirconia incorporated in the model during simulations. The table has been adapted from a previous publication [28] and is being reprinted with the permission of MDPI.

Austenite [26]	Martensite [51]	Zirconia [41]	Units
C11 = 175.0	C11 = 242.3	C11 = 191.0	GPa
C12 = 115.0	C12 = 117.7	C12 = 80.0	GPa
C44 = 135.0	C13 = 45.0	C44 = 40.0	GPa
	C33 = 315.0		GPa
	C44 = 40.5		GPa

.

A simple ductile damage model was incorporated for the austenitic phase [52]. This isotropic ductile damage model is based on the total accumulated plastic slip at a material point, where the local damage ${\phi }_{\mathrm{l}}$ is given by:

$${\phi }_{\mathrm{l}}=min\left(1,\frac{{ϵ}_{crit}}{{{\displaystyle \sum }}_{a=1}^{n_{ss}}{\gamma }^{\alpha }}\right)$$

(9)

The damage phase field value ${\phi }_{\mathrm{l}}=0$ corresponds to fully degraded material stiffness, whereas ${\phi }_{\mathrm{l}}=1$ refers to the fully coherent bulk material. The value of critical plastic strain ${ϵ}_{crit}=0.75$ for the austenite phase was adopted from reported in-situ test results of the similar material [13]. The other parameters and their values are presented in

Table 3. Physical and fitting damage parameter values used for austenitic phase.

Parameter definition	Property	Value	Units
Interface energy	$g_0$	1.0	${\mathrm{Jm}}^{-2}$
Characteristic length	$l_0$	1.0	$\mathsf{\mu }\mathrm{m}$
Damage mobility	$M$	0.01	${\mathrm{s}}^{-1}$
Damage diffusion	$D$	1.0	-
Critical plastic strain	${ϵ}_{crit}$	0.75	-
Damage rate sensitivity	$P$	35	-

. With evolving damage, the element’s stiffness degrades, resulting in no further evolution of physical properties in that element except strain. This false strain in the phase field is not filtered out in the current work as it is negligible compared with the total strain in the material.

2.2.2. An Elastic Model with Brittle Damage for Zirconia Particles

The strength of zirconia particles is many more than the austenite matrix, and therefore for simplicity and computational cost reduction, Zirconia particles are assigned elastic properties only. The elastic stiffness matrix coefficients for the cubic zirconia particles are adopted from the literature. They are shown in Table 2.

A simple brittle damage model was adopted for the zirconia particles. This novel phase-field method (PFM) was derived by Wang et al. [24] for brittle fracture with the crystal plasticity constitutive model. This isotropic brittle model is based on rupture criterion and uses elastic strain energy ${\mathit{W}}_e$ as the driving force. The governing equation of damage evolution is:

$$\mu {\dot{\phi }}_{nl}={\phi }_l-{\phi }_{nl}$$

(10)

where $\mu$ is damage mobility, ${\phi }_l$ is local damage while ${\phi }_{nl}$ is its nonlocal counterpart. The driving force for damage in this model is the elastic strain energy, ${\mathit{W}}_e$, which is calculated as:

$${\mathit{W}}_e=\mathit{E} : {\mathit{S}}^{\ast }$$

(11)

where, $\mathit{E}=\frac{1}{2}\left({\mathit{F}}_e^T{\mathit{F}}_e-\mathit{I}\right)$ is the Green-Lagrange strain and ${\mathit{S}}^{\ast }=ℂ :\mathit{E}$, is its work conjugate, second Piola-Kirchhoff stress. The local damage ${\phi }_l$ is given by:

$${\phi }_l=\frac{{\mathit{W}}_{cr}}{{\mathit{W}}_e}$$

(12)

where,${\mathit{W}}_{cr}$ is the critical strain energy at which damage starts to nucleate. The critical strain energy of $1\times {10}^8 J$ was defined as a damage criterion calculated based on the data available in the literature [53].

2.3. Numerical Simulation and Post-Processing

The simulations are carried out on DAMASK, which provides coupled crystal plasticity and a phase-field damage model. A fast Fourier transformation (FFT) based spectral solver is used to solve mechanical boundary value problems to predict the evolving strain and damage in each phase. Uniaxial monotonic tensile load conditions were defined using the following set of equations:

$${\dot{\mathit{F}}}_{ij}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \ast & 0\\ 0& 0& \ast \end{array}\right]\times {10}^{-3}{ \mathrm{s}}^{-1}$$

(13)

$${\mathit{P}}_{ij}=\left[\begin{array}{ccc}\ast & \ast & \ast \\ \ast & 0& \ast \\ \ast & \ast & 0\end{array}\right] \mathrm{Pa}$$

(14)

where, ${\dot{\mathit{F}}}_{ij}$ and ${\mathit{P}}_{ij}$ are the macroscopic rate of the deformation gradient and macroscopic first Piola–Kirchhoff stress tensor, respectively. The coefficients denoted by ‘*’ highlight the complementary boundary conditions. Using these conditions, the quasi-static uniaxial tensile load was applied in the x-direction, with $1\times {10}^{-3} s^{-1}$ strain rate.

A total strain of 17% was defined in ~850 small increments. Of this, 25 iterations were allowed, with a maximum of 10 cutbacks if needed to find a converged solution for a specific strain increment. Data was recorded for every 5th increment and stored in a binary spectralOut file. With this added complexity, the simulation took ~4000 minutes on a 40 core 100 GB RAM computing cluster with active parallel processing, which is ~650% more computation cost when compared with a similar model without damage definition.

The recorded simulation data was post-processed for obtaining average global and point-based local data for the desired frames. Python scripts and Paraview were extensively used for post-processing the collected data and its visualization in pictorial and graphical form.

3. Results

In this research, a single heterogeneous RVE is considered, and it is used as a case study to run a detailed simulation. The results are processed, broken down into components, and a postmortem analysis is carried out to understand a necessary aspect of the problem at hand. The analysis of the problem is carried out so that the mutual interdependence and/or independence of the underlying physical phenomena can be understood and conclusive relationships can be drawn. This study is aimed to build upon the already available knowledge and information about the deformation behavior of such materials within published literature. It will help make informed decisions about the microstructural engineering of such multi-phase TRIP/TWIP steels in the future.

3.1. Average Global Results

To show simulation results in a simplified manner, each phase’s stress, strain, and damage behavior have been separately processed and averaged for every time increment. The true results for each phase in a comparative manner are shown in Figure 2a. Due to the assignment of elastic properties to zirconia particles, it is observed that brittle damage initiates and propagates aggressively in them after reaching a stress level of 2500 MPa at 2% of global strain. With the propagation of the damage in the zirconia particles, stress relaxation occurs. This brittle damage in the zirconia particles at 2.5% strain or 400 MPa stress has been reported previously by other researchers in the past [13]. For the austenitic phase, it is observed that the yielding occurs at ~200 MPa. After that, a reasonably flat hardening curve follows with ductile damage initiation around 10% of the true global strain. This damage propagation slowly continues to grow with increasing applied external strain.

One might confuse the true behavior of each phase in Figure 2a to presume that the damage in zirconia particles initiates at 2% of the overall strain in the composite, which would not be correct, as to reach this much strain in zirconia particles, the composite material has to undergo higher strain. To clarify this confusion, the normalized stress-strain and damage response of each phase with respect to the composite behavior is provided in Figure 2b. The combined behavior of composite is shown with the help of a solid line in Figure 2, which demonstrates the deformation behavior close to that of the austenitic matrix. This is a reasonable mechanical response as the zirconia particles comprise only 10% of the material composition.

The physics-based TRIP/TWIP model used to simulate the deformation and transformation in the austenitic matrix calculates overall strain accumulation due to different phenomena, which can be disintegrated into its components for study. The breakdown of the overall global material behavior is shown in Figure 3. The evolution of phase transformation, total dislocation density, stress, and strain in the material, including slipping planes with evolving damage, is shown in Figure 3.

With this visualization of property disintegration, the evolution of each attribute and its interdependence on other factors can be noticed. It is observed that the phase transformation initiated at ~4% and exponentially increased up to 12% of true strain, after which the slope of increase becomes linear. Similarly, an increase in the total dislocation density of the material was observed up to ~8% of true strain, after which the slope changes linearly on a logarithmic scale. This drastic increase in the transformation and dislocation density, in the beginning, is due to the large availability of mean free path, which slowly saturates with increasing strain and hence restricts the phase transformation and dislocation motion.

Looking at the material’s stress-strain curve, it is observed that significant brittle damage in the zirconia particles is marked with a kink in the curve, which starts as early as 6.5% of true strain. Despite very early damage initiation and continued propagation, the material continuously endured the applied strain and showed consistent strain hardening up to 13% of true strain. With significant damage evolution in the material, the stiffness degradation and resulting stress relaxation start to appear, which is commonly attributed to necking and material damage. From these results, it is interesting to note that even after the early initiation of damage in the material, phase transformation, dislocation motion, and slight twinning and slipping of planes continues in the rest of the matrix.

This representation of the overall material behavior in Figure 3 depicts the general trends with the help of averaged data. However, the visualization of results at a local scale would help develop a better understanding of why and how these trends occur.

3.2. Local Evolution Maps

The damage initiation and evolution at the local scale are presented in Figure 4a with corresponding stress and strain states. The damage phase field value ${\phi }_l=0$ corresponds to fully degraded material stiffness, whereas ${\phi }_l=1$ refers to the fully coherent bulk material, represented in subsequent figures with black and transparent, respectively.

It is observed that at 9.6% strain and 400 MPa stress, brittle damage in the weakest link of the ceramic particles and ductile damage at the ceramic-matrix interface initiates. These two competing low energy damage phenomena have been reported [13] to initiate at such low strains. With further increase in the strain, up to 13.3%, as shown in Figure 4b, it is observed that the ductile damage on the interface grows 45 degrees to the loading direction, and more brittle damage in zirconia particles evolves. The damage growth and void nucleation lead to the formation of stress and strain localization zones, facilitating further localized damage progression.

At 450 MPa stress and 16.2% overall strain in Figure 4c, it is observed that most of the zirconia particles crack, the void coalescence starts to occur, and the damage progression speeds up. At these large strain fields, the voids also start to initiate in earlier non-damaged austenitic matrix localities. Due to this excessive stiffness degradation of the ceramic and matrix, the load-bearing capacity of the material starts to drop, and the material fails.

The last frame of the converged simulation results is shown in Figure 5 with overlaid damage evolution plane constructions. This figure is drawn according to the deformation gradient vector to closely resemble the actual deformation at each material point.

Ductile damage planes are shown with dotted magenta lines, whereas brittle damage planes are shown with dotted green lines. It is observed that the ductile damage in the austenitic matrix follows the maximum shear plane oriented 45 degrees to the loading axis. On the other hand, the brittle damage in most zirconia particles is observed to have occurred perpendicular to the loading axis. It is also observed that the planes of brittle cracking in zirconia particles intersect the planes of ductile damage progression in the austenitic matrix.

Local distributions of true stress, true strain, martensitic transformation, dislocation density, and triaxial stress are shown in Figure 6. Due to high stiffness, zirconia particles undergo very high stress and low strain. To keep the scales appropriate, the results in Figure 6 are shown only for the steel matrix. Significant heterogeneity in the distribution of these attributes is observed.

Due to high localized stress around the ceramic particles and ceramic/matrix interface, depending on the crystal orientation of austenitic grains, glide systems of higher-order are activated. Hence at lower global strains, enhanced martensite formation is observed in these regions. This martensite transformation of the austenitic phase close to the ceramic particles and the ceramic/matrix interface intensifies with increased global strain, reducing the available mean free path for the dislocation glide; hence the dislocation pinning occurs near these martensitic islands.

This increasing dislocation-pinning results in higher dislocation densities, which further reduces the plastic flow of the material in these zones and results in the local material hardening. These high stress and high strain sites act as potential damage initiation locales where voids start to form and coalesce. While the damage struggles to propagate in these strain-hardened zones, the activation of higher energy slip systems starts taking place in previously less deformed zones, and overall strain in the material continues to grow.

It is observed that the zones with high martensitic volume percentage and high dislocation density also bear the highest amount of stress during deformation, which is in accordance with the previous observations [18,19,20,48]. As the local damage grows in the matrix and ceramic particles, the local stress relaxation occurs due to the loss of stiffness in the damaged elements. It is interesting to note how these phenomena are interconnected at a local scale, shifting from one to the other depending on the amount of activation energy required for a respective phenomenon to take place.

3.3. Local Trends

Local frames are presented in Figure 6, they provide an interesting outlook into the local distribution of attributes. The evolution of these attributes for two different zones of similar size has been recorded and shown in Figure 7. In an interlocked area with expected high strain localization (solid line), the martensitic transformation in the austenite matrix initiates at minor strains. The slope of increasing dislocation density for this region follows a specific trend. Due to this intense increase in the local attributes, damage in this area initiates as early as 10% of the global strain value.

After the excessive strain hardening in these interlocked areas has occurred due to the reduction in the available mean free path, the stress starts to evolve in open areas with relatively low strain localization (dotted outline). The martensitic transformation starts to evolve later but follows the same trajectory. The evolved dislocation density in this area is much lower, with a difference of ~15% at 18% of true strain. This delayed evolution of physical attributes also affects the local damage initiation and evolution scheme.

The line maps of local stress, strain, transformation, and dislocation density at the last frame in highly interlocked and relatively free zones are shown in Figure 8 for a more explicit comparison. For both phases, single-point vertical black lines represent the grain boundaries. It is evident that depending on the grain orientation, the magnitude of attributes for each grain is different. The strain profile in each grain corresponds to the dislocation density. The stress distribution corresponds to the martensitic transformation profile. By visual comparison, it is observed that the local magnitudes of stress, strain, dislocation density, and transformed martensite volume percentage in (i) are significantly lower than in (ii).

Closely looking at the trends in (ii) damage zones and subsequent strain spikes are observed on the ceramic/matrix interface and grain boundaries of highly misoriented grains. The dislocation density, transformation volume percent and subsequent stress increase significantly when moving towards the ceramic/matrix interface. This very high localization of attributes has also been reported by researchers previously [54,55].

4. Discussion

In this work, a physics-based crystal plasticity model incorporated with ductile damage criteria for the meta-stable austenitic matrix and brittle damage criteria for zirconia particles was developed and demonstrated. The model was used to run a full phase simulation on a virtually constructed RVE. Finally, the simulation results were post-processed to analyze the local and global evolution of interdependent deformation and damage phenomena during the tensile deformation of the defined geometry.

The work was carried out by adopting a relatively simplistic approach by assuming the perfect interface between phases. The zirconia particles were only assigned elastic attributes. The geometry was 2D with ~500 grains in total. As indicated in previous publications [37,42], these assumptions resulted in slightly higher contrast in local distributions of the attributes and slightly earlier damage initiation. However, the results are qualitatively accurate and help in establishing the trends which are of interest in the current work.

Weidner et al. [13] performed in-situ tensile deformation of similar material samples to analyze the material damage behavior. The selected results from that publication are shown in Figure 9. The simulation results shown in Figure 5 match well with these experimental observations. Ductile damage in the austenitic matrix is observed to propagate at inclined planes. Brittle damage in the zirconia particles is observed to propagate at perpendicular planes to the loading axis. The stress and strain observed for the damage initiation and propagation also match previously published experimental work [54]. This comparison helps in establishing the validity of the developed model.

From the simulation model results, it is analyzed that there is no evident change observed in the global stress, transformation, dislocation density evolution with the initiation of local damage. The difference in global results is observed at large strain fields, i.e., >20% true strain. Other researchers have also reported this trend in the past [24,31,56]. This is important in understanding that the intended microstructural engineering of such complex multi-phase materials is almost impossible to carry out by postmortem analysis of the failed samples. Instead, the in-situ tests or the phase field models such as the one developed in the current work can help fine-tune the microstructural attributes and verify the intended outcome.

The outlook of this work is to validate the model results accurately by carrying out in-situ tensile tests on well-prepared samples and compare the simulation results with the same area of interest. Another task will be to adopt this model for 3D full phase simulations with varying microstructural attributes for analyzing their effect on mechanical response and development of material formability and process windows.

5. Conclusions

In this work, local damage evolution on micro-mechanical response in TRIP Steel matrix and Mg-PSZ particle composite was analyzed. A combined physical crystal-plasticity-based numerical simulation model with ductile damage for matrix and brittle damage for ceramic particles was developed. Quasi-static monotonic tensile load on a virtual polycrystalline 2D RVE was applied to evaluate the global and local deformation trends and understand the interplay of physical attributes in different regions. From the analysis of the obtained results, the work can be concluded as follows:

• A numerical simulation model based on physics and crystal plasticity for cast X8CrMnNi16-6-6 TRIP steel with ductile damage and elastic cubic zirconia with brittle damage criteria has been successfully developed. It captures the deformation and damage trends in the composite material. Furthermore, the model predicts the necessary artifacts which were observed by researchers experimentally.
• It is observed that with increasing strain, transformation starts at ~5% true strain, and dislocation density saturates at ~14% strain specifically for this material composition. Thus, the damage in different phases initiates at different global strain regimes and affects the localization of the attributes.
• The energy absorbed by the zirconia particles during deformation is relatively high, and significant stiffness degradation in the material is observed due to the damage in this phase. Therefore, the material behavior of composites would be significantly influenced by the percentage of ceramic particles, as their mechanical response will majorly dictate the overall deformation behavior of the composite.
• There is substantial deformation heterogeneity in the matrix depending on the local neighboring conditions. It is observed that at >15% global strain, the interlocked matrix region experiences more significant damage, more stress relaxation, higher dislocation density, reduced mean free path and significantly higher phase transformation than a relatively free matrix region.
• Different grains depict different deformation behavior based on the locality, orientation, and neighboring grains.
• The local damage initiates at the weak junctions of the zirconia particles and the ceramic/matrix interface at as low a global strain as 9%. This local damage keeps evolving without significantly affecting the global stress-bearing capacity of the material. The global stress relaxation starts to occur only after 12% strain.
• The damage in ductile austenite propagates at 45 degrees to the loading axis, and in the zirconia particles, it propagates perpendicular to the loading axis.
• The simulation results are only qualitatively accurate and need further detailed analysis for establishing quantitative accuracy.