Skip to main content
  • Research article
  • Open access
  • Published:

Particle interspacing effects on the mechanical behavior of a Fe–TiB2 metal matrix composite using FFT-based mesoscopic field dislocation mechanics

Abstract

This paper presents an application to metal matrix composites (MMCs) of an enhanced elasto-viscoplastic Fast Fourier Transform (EVP-FFT) formulation coupled with a phenomenological continuum Mesoscale Field Dislocation Mechanics (MFDM) theory. Contrary to conventional crystal plasticity, which only accounts for plastic flow and hardening induced by statistically stored dislocations (SSDs), MFDM-EVP-FFT also describes the evolution of polarized geometrically necessary dislocation (GND) density and its effect on both plastic flow and hardening. Numerical results for a Fe–TiB2 MMC made of a ferrite matrix (α-Fe) and elastic ceramic particles (TiB2) are presented. Full-field simulations are performed using synthetic periodic unit cells representative of the MMC, with single-crystalline and polycrystalline matrix, for different particle interspacing distances. A strong dependence of the predicted equivalent stress, cumulated plastic strain and GND density fields with particle interspacing distance is observed, in contrast with conventional crystal plasticity. Correlations between these mechanical fields and microstructural features, and their influence on local and global mechanical behavior are examined for the different MMC microstructures.

Introduction

The influence of the size of microstructural features on the macroscopic mechanical behavior of crystalline materials has been the subject of many studies over past decades. These size effects are controlled by the development of plastic strain gradients in the microstructure. A decrease in characteristic sizes of the microstructure may lead to an increase of the overall stress for a given applied strain, and therefore contribute to the material’s hardening. This tendency is often described as a smaller is stronger size effect.

For metallic polycrystalline materials, grain-size strengthening is one the most common mechanisms. This size effect is also known as the Hall–Petch effect [1, 2], which states that the flow stress at constant strain is proportional to the inverse square root of the mean grain size. This grain-size dependence is related to the development of plastic strain gradients, mainly due to lattice incompatibilities associated to the mismatch of plastic deformation between neighboring grains. The concept of geometrically necessary dislocations (GNDs) was introduced to rationalize the accommodation of plastic deformation gradients in these regions [3]. GNDs were introduced to capture size effects, and were distinguished from the statistically stored dislocations (SSDs) that are trapped within grains and contribute to size-independent hardening only. As the global strain level increases, strain gradients are more pronounced in the vicinity of grain boundaries where GND pile-ups are preferentially formed. Hence, polycrystals with smaller grain sizes—and therefore with higher GND content—exhibit a stronger mechanical response. This behavior excludes the case of nano-crystalline materials with sufficiently small grain sizes, for which a reverse Hall–Petch effect was observed, see e.g. [4].

Similarly, the macroscopic behavior of polycrystalline materials can be influenced by the morphology of grains, and, in the case of multi-phase aggregates, by the volume fraction and the spatial distribution of the different phases. Inhomogeneous deformation processes, such as the development of highly deformed slip bands, also have a strong influence on the macroscopic response of polycrystalline aggregates. In particle-reinforced metals, inhomogeneous deformation involving strain accommodation between the soft matrix and the hard particles is at the origin of strain gradients. Therefore, in such materials, additional characteristic length-scales can be associated with the size of the reinforcements and the distance between particles. In general, flow stress of composites increases with decreasing particle size at constant volume fraction of reinforcements [5, 6].

In conventional crystal plasticity formulations, customarily used to model single-crystal and polycrystal behavior, intrinsic material length-scales are not considered. Strain hardening only depends on crystallographic slip variables, and effects of GNDs are ignored. Therefore, these models are unable to capture size effects. Theories based on strain gradient plasticity have been developed in the last 3 decades to overcome this limitation [7,8,9,10,11,12,13,14]. Many models directly incorporate GND densities in the evolution of the flow stress leading to reported extended hardening laws containing both GND and SSD densities, see e.g. [7, 15,16,17,18].

In the framework of the elastic theory of continuously distributed dislocations [19], Acharya and co-workers [20,21,22] developed a continuum dislocation mechanics theory, called Field Dislocation Mechanics (FDM), based on the use of Nye’s dislocation density tensor [23] as an internal variable. FDM enables prediction of the time-dependent mechanical response of crystalline solids at scales in which the distribution of dislocations can be represented by GND densities. Subsequently, an approach based on space–time averaging of the FDM field equations, called Phenomenological Mesoscopic Field Dislocation Mechanics (PMFDM) was proposed [24]. PMFDM extends the theory at a mesoscopic scale, integrating the mobilities of both SSDs and GNDs and their contribution to strain hardening through a phenomenological specification of the velocity of polarized dislocations.

A finite element (FE) based implementation of PMFDM was first presented in [25, 26]. Since then, several FE numerical studies using PMFDM have been performed to analyze size effects in single crystals [27, 28] and multicrystalline thin films [29, 30], including the role of GNDs on the directionality of the yield stress in strain-aged steels [31]. Furthermore, a numerical implementation of a reduced version of the PMFDM theory was presented in [32] and used to study the work hardening of Al-based metal-matrix composites in [33]. In general, PMFDM results were consistent and in qualitative agreement with experimental observations.

As an efficient alternative to FE for numerically solving micromechanical problems, Fast Fourier Transform (FFT)-based methods were initially developed and applied to composites [34,35,36,37], in which the material’s heterogeneity is given by the spatial distribution of phases with different mechanical properties. FFT-based methods were later adapted to polycrystalline materials [38,39,40], where the heterogeneity is related to the spatial distribution of anisotropic crystals with different orientations. The original crystal plasticity FFT-based implementations showed the feasibility of efficiently solving the micromechanical behavior of polycrystalline unit cells with complex microstructures. In the context of continuum dislocation mechanics, the FFT-based method was adapted and applied to solve the elastostatic field equations of the FDM theory [41,42,43]. In addition, in the dynamic regime, the GND density transport equation was solved using a spectral exponential filter to avoid numerical oscillations related to Gibbs phenomenon [44]. Furthermore, the conventional (based on local crystal plasticity, no GND transport, etc.) elasto-viscoplastic FFT formulation [39], denoted CP-EVP-FFT here, was adapted to consider PMFDM equations in its reduced version. This implementation, abbreviated MFDM-EVP-FFT, was originally developed and demonstrated for two-phase laminate composites by [45] and applied to Voronoi polycrystals with different average grain sizes under monotonic tensile and reversible tension–compression loadings to study the Bauschinger effect [46]. Compared to FE, FFT-based methods may be advantageous to avoid mesh generation [47] and less time-consuming (in some major cases except exascale calculations) [48].

This contribution presents a new application of MFDM-EVP-FFT to study the mechanical behavior of a new type of metal matrix composites (MMC), made of a full ferritic matrix (α-Fe, body-centered cubic (BCC) crystallographic structure) and ceramic particles of type TiB2 (hexagonal compact (HCP) crystallographic structure). This new composite is attractive for industrial applications due to the very high elastic moduli and low density of TiB2 particles, leading to an overall increase in the specific stiffness of at least 15% [49]. Previous works investigated the structure and chemistry of interfaces [50] and the damage mechanisms of Fe–TiB2 [51, 52]. It was demonstrated that cracking of TiB2 particles and subsequent damage initiation in the adjacent matrix occur, resulting in a significant loss of effective moduli, toughness and ductility, which is detrimental to industrial applications of such composites. In addition, the partitioning of stresses and strains between the ductile ferrite matrix and the hard ceramic reinforcements was shown to depend on particle sizes, shapes, distribution and volume fraction [53]. Therefore, these microstructural parameters are likely to alter the deformation characteristics in a significant way, and therefore their influence on the mechanical behavior needs to be assessed in order to optimize the material’s processing routes.

The objective of this paper is to use the MFDM-EVP-FFT formulation in order to describe internal size effects on the mechanical behavior of the Fe–TiB2 composite. Specifically, particle interspacing effects and their influence on the development of stresses, plastic distortion and distribution of dislocation densities will be discussed. For this purpose, 3D synthetic microstructures are used as input of FFT-based simulations.

The paper is organized as follows. In “Description of MFDM-EVP-FFT model and input material data” section, the constitutive equations of the MFDM and the elasto-viscoplastic FFT-based numerical implementation for MFDM are presented, as well as the material’s microstructure and parameters used as input for numerical simulations. “Effect of the particle interspacing distance” section shows the application of the MFDM-EVP-FFT model to synthetic microstructures with single-crystal matrix and different interspacing distances between TiB2 particles to study the influence of particle interspacing on mechanical behavior of these composites. In “Interplay between particle interspacing and grain size effects” section, a polycrystalline ferritic matrix is considered to study the interplay between grain size and particle interspacing. In “Conclusions” section we present our conclusions.

Description of MFDM-EVP-FFT model and input material data

MFDM equations

To solve the displacement \( \varvec{u} \) and stress \( \varvec{\sigma} \) fields, the following equations are solved in a small deformation setting for elasto-viscoplastic behavior in a volume \( V \) with external boundary \( S \) using standard traction/displacement boundary conditions on \( S_{t} \) and \( S_{u} \) (\( S = S_{t} \cup S_{u} \)):

$$ \begin{aligned} {\text{div}} \varvec{\sigma}= 0\; \hfill \\\varvec{\sigma}= {\mathbf{C}}:\left( {{\mathbf{U}}^{{\mathbf{e}}} } \right)^{sym} \hfill \\ {\mathbf{U}} = {\text{grad }}{\mathbf{u}} = {\mathbf{U}}^{{\mathbf{e}}} + {\mathbf{U}}^{{\mathbf{p}}} \hfill \\\varvec{\sigma}\cdot {\mathbf{n}} = \overline{{\mathbf{T}}}\; {\text{on}} \;S_{t} \hfill \\ {\mathbf{u}} = \overline{{\mathbf{u}}}\; {\text{on}}\,S_{u} \hfill \\ \end{aligned} $$
(1)

where \( {\mathbf{C}} \) is the fourth order elastic stiffness tensor with classic minor and major symmetries.

In the presence of dislocation ensembles [20, 24,25,26], both the average plastic distortion \( {\mathbf{U}}^{{\mathbf{p}}} \), which results from dislocation motion, and the average elastic (or lattice) distortion \( {\mathbf{U}}^{{\mathbf{e}}} \) are incompatible fields. In non-local crystal plasticity, and depending on the resolution scale, dislocation ensembles can be categorized as GND and SSD [3]. Indeed, if the resolution scale is fine enough (microscopic scale), all dislocations appear to be polarized, i.e. GNDs. At a higher resolution (mesoscopic scale), both SSDs and GNDs are distinguished. SSDs accumulate inside grains and only contribute to plastic flow, whereas GNDs contribute to both plastic flow and long-range internal stresses. Polarized dislocations are represented by the dislocation density (or Nye) tensor \( \varvec{\alpha} \). The mesoscale FDM theory is based on an averaged value of the tensor \( \varvec{\alpha}. \) Here, a simplified reduced version of the MFDM is considered [24, 32], where incompatible fields are assumed to be as smooth as necessary, and the average plastic distortion rate writes:

$$ {\dot{\mathbf{U}}}^{{\mathbf{p}}} =\varvec{\alpha}\times {\mathbf{v}} + {\mathbf{L}}^{{\mathbf{p}}} . $$
(2)

The mobility of SSDs is represented by the mesoscale plastic distortion rate \( {\mathbf{L}}^{{\mathbf{p}}} \), where the averaging procedure was defined in [24]. The space–time evolution of the average dislocation density tensor \( \varvec{\alpha} \) is obtained from the conservation of dislocation flux and is prescribed as [54, 55]:

$$ \dot{\varvec{\alpha }} = - {\text{curl}}\left( {{\dot{\mathbf{U}}}^{{\mathbf{p}}} } \right). $$
(3)

Constitutive specifications on the dislocation velocity \( {\mathbf{v}} \) and on the slip distortion rate \( {\mathbf{L}}^{{\mathbf{p}}} \) are given from thermodynamic considerations following the theory introduced by Acharya and Roy [24]. Furthermore, plastic flow incompressibility is considered. The GND velocity \( {\mathbf{v}} \) is prescribed as follows:

$$ {\mathbf{v}} = \frac{{\mathbf{g}}}{{\left| {\mathbf{g}} \right|}}v\quad{\text{with}}\;v \ge 0 $$
(4)

where \( {\mathbf{g}} \) is the glide force parallel to \( {\mathbf{v}} \) and \( v \) is the magnitude of \( {\mathbf{v}} \). The constitutive equation adopted for \( v \) is based on the Orowan law for GND mobile dislocations together with a Taylor hardening rule. Here, a mechanistic formulation for \( v \) similar to [29] is used:

$$ v = \frac{{\zeta^{2} b}}{N}\left( {\frac{\mu }{{{{\uptau }}_{c} }}} \right)^{2} \mathop \sum \limits_{{{\text{s}} = 1}}^{N} \left| {{\dot{\upgamma}}^{\text{s}} } \right| \ge 0 $$
(5)

where \( N \) is the total number of slip systems (\( N = 24 \) for the BCC structure in ferrite: α-Fe), \( {\dot{\upgamma}}^{\text{s}} \) is the slip rate on slip system \( s \), \( \upzeta \) is a material constant, \( b \) is the magnitude of the Burgers vector, \( {\uptau}_{c} \) is the shear strength and \( \mu \) is the isotropic elastic shear modulus of the material.

In crystal plasticity, the plastic distortion rate tensor \( {\mathbf{L}}^{{\mathbf{p}}} \) due to slip is defined as:

$$ {\mathbf{L}}^{{\mathbf{p}}} = \mathop \sum \limits_{{{\text{s}} = 1}}^{N} {\dot{\upgamma}}^{\text{s}} {\mathbf{b}}^{{\mathbf{s}}} \otimes {\mathbf{n}}^{{\mathbf{s}}} $$
(6)

where \( {\mathbf{b}}^{{\mathbf{s}}} \) and \( {\mathbf{n}}^{{\mathbf{s}}} \) denote, respectively for each slip system \( s \), the slip direction and the slip plane unit normal. The constitutive equation for \( {\dot{\upgamma}}^{\text{s}} \) introduced in Eqs. 5 and 6 is given by a classic viscoplastic flow rule as a power law:

$$ {\dot{\upgamma}}^{\text{s}} = {\dot{\upgamma}}^{0} \left( {\frac{{\left| {{\uptau}^{\text{s}} } \right|}}{{{\uptau}_{\text{c}} }}} \right)^{{1/{\text{m}}}} {\text{sign}}\left( {{\uptau}^{\text{s}} } \right) $$
(7)

where \( {\text{m}} \) is the strain rate sensitivity of the material, \( {\uptau}^{\text{s}} \) is the resolved shear stress, \( {\dot{\upgamma}}^{0} \) is the reference slip rate and \( {\uptau}_{\text{c}} \) is considered identical for all slip systems. As numerical applications in “Interplay between particle interspacing and grain size effects” section are not concerned by cyclic plasticity, no intra-crystalline phenomenological back-stress evolution has been introduced (as it was in [29]). The cumulated slip rate on all slip systems due to both GNDs and SSDs is given by:

$$ {\dot{{\Gamma }}} = \left| {\varvec{\alpha}\times {\mathbf{v}}} \right| + \mathop \sum \limits_{{{\text{s}} = 1}}^{N} \left| {{\dot{\upgamma}}^{\text{s}} } \right|. $$
(8)

The evolution law for the shear strength \( {\uptau}_{\text{c}} \) follows the same hypotheses as the strain-hardening model developed by [29], which is an extension of earlier models derived by [17, 56]:

$$ {\dot{\uptau}}_{c} = {{\uptheta}_{0}} \frac{{{\uptau}_{\text{s}} - {\uptau}_{\text{c}} }}{{{\uptau}_{\text{s}} - {\uptau}_{0} }}{\dot{{\Gamma }}} + k_{0} \frac{{{\zeta}^{2} \mu^{2} b}}{{2\left( {{\uptau}_{\text{s}} - {\uptau}_{0} } \right)}}\left( {\mathop \sum \limits_{{{\text{s}} = 1}}^{N} \left| {\varvec{\alpha}\cdot {\mathbf{n}}^{{\mathbf{s}}} } \right|\left| {{\dot{\upgamma}}^{\text{s}} } \right| + \mathop \sum \limits_{{{\text{s}} = 1}}^{N} \left| {\varvec{\alpha}\cdot {\mathbf{n}}^{{\mathbf{s}}} } \right|\left| {\varvec{\alpha}\times {\mathbf{v}}} \right|} \right) $$
(9)

where \( \uptau_{0} \) is the yield strength due to lattice friction, \( {\uptau}_{\text{s}} \) is the saturation stress, \( \uptheta_{0} \) is the initial hardening rate and \( k_{0} \) is a parameter related to a geometric mean free path due to GND forest on slip system s [17]. In the case of a model based on conventional crystal plasticity (no GND, i.e. \( \varvec{\alpha}= 0 \)), Eq. 9 reduces to a classic Voce-Kocks law [57].

Numerical implementation of MFDM-EVP-FFT

In the FFT-based framework for periodic heterogeneous media, the constitutive equations in (1) can be solved using the Green’s function method [58] through an integral Lippmann–Schwinger equation. The Fourier transform of this equation is expressed as follows:

$$ \begin{aligned} \widehat{\varvec{\varepsilon}}\left(\varvec{\xi}\right) &= - \widehat{{\varvec{\Gamma}}}^{0} \left(\varvec{\xi}\right):\widehat{\varvec{\tau}}\left(\varvec{\xi}\right) \varvec{ }\forall\varvec{\xi}\ne 0 \\ \widehat{\varvec{\varepsilon}}\left( 0 \right) & = \varvec{ \varepsilon } \end{aligned}$$
(10)

where \( \varvec{\xi} \) is the Fourier vector, \( {\varvec{\Gamma}}^{0} \) is the modified Green tensor associated with the homogeneous elastic moduli \( {\mathbf{C}}^{0} \), and \( \tau_{ij} = \sigma_{ij} - C_{ijkl}^{0} \varepsilon_{kl} \) represents the stress polarization tensor field due to heterogeneities. As for notation, \( \widehat{\varvec{\varepsilon}}\left(\varvec{\xi}\right) \) is the continuous Fourier transform of tensor field \( \varvec{\varepsilon}\left( \varvec{x} \right) \). The calculation of the modified Green tensor in the Lippmann–Schwinger equation is performed using a centered finite difference scheme on a rotated grid introduced in [59].

Equation 10 is solved using an augmented Lagrangian iterative scheme introduced in [37]. Stress and strain fields at iteration \( \left( {\text{n}} \right) \) are approximated by auxiliary fields \( \lambda_{ij}^{{\left( {\text{n}} \right)}} \) and \( e_{kl}^{{\left( {\text{n}} \right)}} \) respectively. The stress polarization tensor then becomes:

$$ \tau_{ij}^{{\left( {\text{n}} \right)}} = \lambda_{ij}^{{\left( {\text{n}} \right)}} - C_{ijkl}^{0} e_{kl}^{{\left( {\text{n}} \right)}} . $$
(11)

The augmented Lagrangian scheme also requires the nullification of the residual \( \varvec{R} \), which depends on the stress and strain tensors \( \varvec{\sigma}^{{\left( {{\text{n}} + 1} \right)}} \) and \( \varvec{\varepsilon}^{{\left( {{\text{n}} + 1} \right)}} \) (Eq. 12). This non-linear equation is solved using a Newton–Raphson iterative scheme.

$$ R_{ij} \left( {\varvec{\sigma}^{{\left( {{\text{n}} + 1} \right)}} } \right) = \sigma_{ij}^{{\left( {{\text{n}} + 1} \right)}} + C_{ijmn}^{0} \varepsilon_{mn}^{{\left( {{\text{n}} + 1} \right)}} \left( {\varvec{\sigma}^{{\left( {{\text{n}} + 1} \right)}} } \right) - \lambda_{ij}^{{\left( {\text{n}} \right)}} - C_{ijmn}^{0} e_{mn}^{{\left( {{\text{n}} + 1} \right)}} . $$
(12)

Once the convergence is reached for \( \varvec{\sigma}^{{\left( {{\text{n}} + 1} \right)}} \) and \( \varvec{\varepsilon}^{{\left( {{\text{n}} + 1} \right)}} \), the new guess for the auxiliary stress field \( \varvec{\lambda} \) is given by using Uzawa’s descent algorithm:

$$ \lambda_{ij}^{{\left( {{\text{n}} + 1} \right)}} = \lambda_{ij}^{{\left( {\text{n}} \right)}} + C_{ijkl}^{0} \left( {e_{kl}^{{\left( {{\text{n}} + 1} \right)}} - \varepsilon_{kl}^{{\left( {{\text{n}} + 1} \right)}} } \right). $$
(13)

The iterative algorithm is stopped when the normalized average differences between the stress fields \( \varvec{\sigma} \) and \( \varvec{\lambda} \) and the strain fields \( \varvec{\varepsilon} \) and \( \varvec{e} \) are smaller than a threshold error of \( 10^{ - 5} \). This condition implies the fulfillment of both stress equilibrium and strain compatibility up to sufficient accuracy.

In the algorithm described above, an overall macroscopic strain \( {\text{E}} = \varepsilon^{{\left( {\text{n}} \right)}} \) is applied to the periodic unit cell \( V \) in the form of:

$$ E_{ij} = E_{ij}^{t} + \dot{E}_{ij} {\Delta}t. $$
(14)

In cases of mixed boundary conditions with imposed macroscopic strain rate \( \dot{E}_{ij} \) and stress \( {\varSigma}_{ij} \), the (n + 1)-guess of the macroscopic strain \( E_{ij}^{{\left( {{\text{n}} + 1} \right)}} \) was given in [37, 39].

For the space–time evolution of the dislocation density tensor (Eq. 3), an explicit forward Euler scheme was derived to numerically solve this equation starting from an implicit backward Euler scheme together with a Taylor expansion at first order of \( \alpha_{ij}^{{t + {{\Delta}}t}} \) where \( {{\Delta }}t \) is the time step. An efficient numerical spectral approach that uses an exponential low-pass filter [44] in order to avoid non accurate and unstable solutions due to the occurrence of high-frequency Gibbs oscillations and due to numerical instabilities resulting from the hyperbolic nature of Eq. 3. The exponential low-pass filter is defined as function of frequencies \( {\upeta} \) as:

$$ {{\upkappa }}\left( {\upeta} \right) = { \exp }\left( { - {{\upbeta }}\left( {\upeta} \right)^{{2{\text{p}}}} } \right) $$
(15)

where \( {{\upbeta }} = - { \log }\left( {{{\varepsilon }}_{M} } \right) \) and \( 2{\text{p}} \) correspond to the high-frequency damping parameter and the order of the filter respectively. Applying the filter to Eq. 3 leads to:

$$ \widehat{\alpha }_{ij}^{{t + {{\Delta }}t}} = {{\upkappa }}\left( {\upeta} \right)\left( {\widehat{\alpha }_{ij}^{t} - {{\Delta }}t {\text{i}} \xi_{k} \left( {\left( {\widehat{{\alpha_{ij} v_{k} }}} \right)^{t} - \left( {\widehat{{\alpha_{ik} v_{j} }}} \right)^{t} } \right)} \right) - {{\Delta }}t {\text{i }}\xi_{k} e_{jkl} \left( {\widehat{{{\text{L}}_{il}^{\text{p}} }}} \right)^{t} . $$
(16)

For the second term in Eq. 15 not concerned by the filter, a centered finite difference scheme coupled with discrete Fourier transform is used [45]. To satisfy stability requirements for the numerical resolution of the dislocation density transport equation, the time step \( {{\Delta }}t \) is fixed according to a Courant–Friedrichs–Lewy (CFL) condition:

$$ {{\Delta }}t_{CFL} = c\frac{\delta }{{v_{max} }} $$
(17)

where \( {{\Delta }}t_{CFL} \) is the refined time step, \( \delta \) is the voxel size, \( c \) is a user-specified constant and \( v_{max} \) is the maximal GND velocity.

Material and simulation data

Description of material

The material studied here is a ferritic-steel-based composite reinforced with a 9.5% volume fraction of TiB2 particles. The composite was produced by ArcelorMittal Research SA, by ingot casting with in situ precipitation of the ceramic TiB2 particles during solidification by eutectic reactions [49]. The products were cold rolled to produce 1.5 mm thick sheet material before being subjected to an annealing treatment. As result of this processing route, the α-ferrite grains have a body-centered cubic structure (BCC) and presents no crystallographic texture. The TiB2 particles are uniformly distributed in the matrix, elongated preferentially in the rolling direction and exhibit hexagonal-like sections. Moreover, the crystal structure of the TiB2 particles is hexagonal—as it was determined by X-ray powder diffraction analysis [49]—with a preferential alignment of the rolling direction with the c-axis of the hexagonal structure. The mean grain size of the matrix is about 4.5 µm and the mean equivalent diameter of TiB2 particles is close to 1.7 µm. These descriptors of shapes and orientations of grains and particles were derived from statistical analysis of 2D SEM images of the material (not documented in the present paper).

Material and simulation parameters

This section presents numerical simulations performed on different periodic unit cells representative of the Fe–TiB2 MMC. These synthetic unit cells are described in “Effect of the particle interspacing distance” and “Interplay between particle interspacing and grain size effects” sections and consist of an elasto-viscoplatic matrix with purely elastic particles. Assumption of elastic behavior for particles is consistent with the strong stiffness contrast with the ferrite matrix, considering the low total strain level reached in the simulations (up to 0.2%).

Due to periodic boundary conditions required by FFT-based methods, periodicity of the simulation volume is assumed in all three spatial directions. Elastic constants for the ferrite matrix and the particles are reported in Table 1 and reflect the cubic and hexagonal structures of α-ferrite and TiB2, respectively. For the elasto-viscoplatic matrix, the material parameters related to slip rule and GND velocity (\( \dot{\gamma }_{0} , m \) in Eq. 7 and \( \zeta \) in Eqs. 5 and 9) are consistent with pure α-Fe and are given in Table 2. A specific fit to experimental data was carried out only to identify material parameters related to the hardening model (\( \tau_{0} , \tau_{s} \) and \( \theta_{0} \) in Eq. 9). The identification procedure is described in “Identification of hardening parameters” section. Material parameter \( k_{0} \) has not been experimentally identified and is fixed following a value chosen by [17, 25], i.e. \( k_{0} = 20 \) in Eq. 9. The Burgers vector magnitude for α-Fe is \( b = 2.48 \times 10^{ - 10} m \).

Table 1 Elastic constants used for numerical simulations [60]
Table 2 Material parameters related to plasticity used for numerical simulations

Given the body-centered cubic structure of the α-Fe lattice, both \( \left\{ {110} \right\} \left\langle {\overline{1} 11} \right\rangle \) and \( \left\{ {211} \right\} \left\langle {\overline{1} 11} \right\rangle \) slip systems are considered (i.e. a total of 24 slip systems, with \( \left\{ {123} \right\} \left\langle {\overline{1} 11} \right\rangle \) excluded in this analysis). As the volume fraction effect is not investigated in the present paper, the volume fraction of particles is fixed at 9.5%, with slight variations, depending on the chosen voxelization. In the following, all considered unit cells are subjected to a pure uniaxial tensile loading in the Z-direction, with mixed strain/stress boundary conditions and applied strain rate \( \dot{E}_{ZZ} = 0.001 \; {\text{s}}^{ - 1} \). No specific condition was considered across material interfaces in this study, though such physically motivated conditions can be introduced in the MFDM theory [29]. Finally, the specific numerical parameters used for GND density transport equation in Eqs. 15 and 17 are taken from a previous study [44], where these parameters were optimized: \( p = 1, {{\varepsilon }}_{M} = 0.2, c = 0.25 \).

Identification of hardening parameters

A fit to experimental tensile data is performed to identify material parameters related to the hardening model (\( \tau_{0} , \tau_{s} \) and \( \theta_{0} \) described in “MFDM equations” section). To this end, a periodic unit cell representative of the Fe–TiB2 composite was generated using DREAM3D software [61] and is shown in Fig. 1a. The unit cell is composed of about 200 equiaxed grains and a 9.5% volume fraction of uniformly distributed particles in the form of hexagonal cylinders. Ferrite grains have a random crystallographic texture and their sizes follow a log-normal distribution with mean grain size of \( \overline{{d_{g} }} = 4.5\;\;\upmu{\text{m}} \) and standard deviation of about 0.4 × \( \overline{{d_{g} }} \). TiB2 particles have a fixed aspect ratio of 3 with the cylinder axis in the Z-direction. Their equivalent diameter also follows a log-normal distribution with mean size of \( \overline{{d_{p} }} = 1.7\;\upmu{\text{m}} \) and standard deviation of about 0.6 × \( \overline{{d_{g} }} \). Finally, the c-axis of the HCP lattice is oriented along Z-axis for all particles.

Fig. 1
figure 1

a RVE of the two-phase composite discretized with 64 × 64 × 64 voxels and used for identification of hardening parameters; b macroscopic stress–strain responses issued from experimental tensile test (dashed black) and obtained from numerical simulation with fitted parameters and RVE of a

Parameters have first been fitted using a coarse 32 × 32 × 32 discretization and then adjusted for the final more refined 64 × 64 × 64 discretization. Several realizations were considered to ensure representative numerical results. Both the experimental tensile stress–strain response and the numerical simulation with the identified hardening parameters are displayed in Fig. 1b.

Outputs

For all performed numerical simulations, overall stress/strain responses are reported and mechanical fields inside of grains are investigated. The latter is based on the analysis of three different scalar outputs recorded at the same overall tensile strain \( E_{ZZ} = 0.2\% \):

  1. i.

    The equivalent Von Mises stress \( \sigma_{eq} \), in order to evidence stress hotspots in the microstructures:

    $$ \sigma_{eq} = \sqrt {\frac{3}{2}{\mathbf{s}} :{\mathbf{s}}} $$
    (18)

    where \( {\mathbf{s}} \) is the deviatoric stress tensor.

  2. ii.

    The equivalent cumulated plastic strain \( \varepsilon_{eq}^{p} \), to analyze plastic strain localization in the microstructures:

    $$ \varepsilon_{eq}^{p} = \mathop \int \limits_{0}^{t} \sqrt {\frac{2}{3}\dot{\varvec{\varepsilon }}^{\varvec{p}} :\dot{\varvec{\varepsilon }}^{\varvec{p}} } dt $$
    (19)

    where \( \dot{\varvec{\varepsilon }}^{\varvec{p}} \) is the symmetric part of \( \dot{\mathbf{U}}^{\mathbf{p}} \).

  3. iii.

    The L2-norm \( {\varvec{\upalpha}} \) of the Nye tensor or alternatively the GND density, defined as:

    $$ \rho_{GND} = \frac{{\varvec{\upalpha}}}{b} = \frac{{\sqrt {{\varvec{\upalpha}}:{\varvec{\upalpha}}} }}{b}. $$
    (20)

The simplified analysis of \( {\varvec{\upalpha}} \) provides a scalar measure of GND, instead of the more involved study of each of the 9 Nye tensor components.

Effect of the particle interspacing distance

In this section we consider a 64 × 64 × 64 two-phase synthetic unit cell composed of one elastic particle embedded in a single crystal elasto-viscoplastic matrix, see Fig. 2a. Motivated by experimental observations, the geometry of the particle is a hexagonal prism, as illustrated in Fig. 2a. Given the periodic boundary conditions, this unit cell corresponds to a periodic cubic distribution of particles in an infinite single crystal matrix. Hence, spacings \( s_{X} \), \( s_{Y} \) and \( s_{Z} \) between particles are constant in each direction but slightly vary from one direction to another (\( s_{X} \ne s_{y} \ne s_{z} \)) owing to the hexagonal prismatic shape of particles. We here limit the study to the effect of the mean particle interspacing—denoted \( l_{p} \)—considered as a single internal length scale parameter and defined as the averaged value:

$$ l_{p} = \frac{1}{3}\left( {s_{X} + s_{Y} + s_{Z} } \right). $$
(21)
Fig. 2
figure 2

a Single-particle two-phase synthetic unit cell discretized with 64 × 64 × 64 voxels; b studied crystallographic orientations of the ferrite matrix displayed in the standard stereographic triangle

The choice of a periodic cubic distribution of particles here ensures that \( l_{p} \) is constant. It should be noted that the aspect ratio of the particle (close to 1) was chosen lower than experimentally observed particles, so that no strong variation between \( s_{X} \), \( s_{Y} \) and \( s_{Z} \) is introduced.

Four different realistic mean particle interspacings ranging from \( l_{p} = 0. 10 5\;\upmu{\text{m}} \) to \( l_{p} = 1.26\;\upmu{\text{m}} \) are studied, as well as two different crystallographic orientations. The latter are considered to provide some qualitative insight of the effect of crystallographic orientation on macroscopic response and local mechanical fields. The first one is associated to multiple slip and such that \( \left[ {\overline{1} 11} \right] \) lies parallel to the tensile direction Z. For the second, single slip is the preferred slip mode (with \( \left( {1\overline{1} 2} \right)\left[ {\overline{1} 11} \right] \) system exhibiting the highest Schmid factor) and \( \left[ {\overline{5} ,6,27} \right] \) is parallel to the tensile direction Z. These orientations will be referred to as orientation #1 and orientation #2, respectively. Finally, the MFDM-EVP-FFT predictions are compared to corresponding results obtained with the conventional CP-EVP-FFT formulation of [39].

Effect of interspacing distance on macroscopic response

Macroscopic strain–stress responses of the synthetic RVEs for the range of investigated mean particle interspacing distances are reported on Fig. 3a, b for both crystallographic orientations. For both orientations, an increase of the flow stress is observed with interspacings decreasing from \( l_{p} = 1.26\;\upmu{\text{m}} \) down to \( l_{p} = 0.105\;\upmu{\text{m}} \). This macroscopic smaller is stronger effect of particle interspacing was already obtained with strain gradient plasticity models applied to MMCs, see e.g. FEM-based calculations of [62] or [63] who used the model of [12, 13]. Higher stress magnitudes are reached for orientation #1 compared to orientation #2, which is consistent with the lower Schmid factors of slip systems in the latter case, i.e. orientation #1 is harder with respect to the applied tensile direction (the maximum Schmid factor is 0.314). The flow stress at 0.2% macroscopic strain corresponding to the lowest \( l_{p} \) is 157% higher than the initial macroscopic yield stress for orientation #1 and 130% higher for orientation #2. Therefore, as it is directly observed from Fig. 3a, b, the macroscopic tensile stress is more sensitive to the particle interspacing effect for the multiple slip orientation #1 than for the single slip orientation #2. No scaling law of the macroscopic flow stress has been identified here, in contrast with grain-size effect for which the well-known Hall–Petch’s law was obtained from other simulations, not reported here, see [46].

Fig. 3
figure 3

Macroscopic stress–strain responses of the two-phase single-particle unit cell under uniaxial tension predicted by the MFDM-EVP-FFT formulation (solid lines) with interspacing distances \( l_{p} \) ranging from 0.105 to 1.26 µm for the two orientations of the single crystal matrix, orientations #1 (a) and #2 (b), respectively. For comparison, tensile responses given by conventional plasticity (CP-EVP-FFT) are also shown (dotted line)

In contrast with the MFDM-EVP-FFT predictions, results obtained with the conventional CP-EVP-FFT model are size independent. Also, the macroscopic hardening is significantly lower in conventional results. This is due to the absence of description of GND densities, that contribute to the overall hardening in the MFDM-EVP-FFT approach, as will be seen in the following paragraph.

Effect of interspacing distance on intra-granular mechanical fields

The fields of scalar quantities defined in “Outputs” section are analyzed here to study stress and strain localization and distribution of dislocation densities. The mechanical fields correspond to the same macroscopic tensile strain \( E_{ZZ} = 0.2\% \).

The spatial distribution of \( \sigma_{eq} \) in different sections of the unit cell are reported in Fig. 4. As expected, higher stresses are found in the elastic particles. As illustration, the average \( \sigma_{eq} \) particle-to-matrix ratio is around 1.7 for orientation #1 and 2.3 for orientation #2 in the case with the largest \( l_{p} \). The difference in stress amplitudes between the two studied crystallographic orientations mentioned in “Effect of interspacing distance on macroscopic response” section is clearly observed in Fig. 4a–d.

Fig. 4
figure 4

ad Spatial distributions of equivalent Von Mises stress \( \sigma_{eq} \) predicted with MFDM-EVP-FFT for a macroscopic strain \( E_{ZZ} = 0.2\% \), for two different particle interspacing distances \( l_{p} \) and two different crystallographic orientations of the matrix. 2-D sections containing the unit cell center are either parallel or perpendicular to the Z tensile direction (used axis systems are displayed in a). Each 2-D section represents the same physical length. For comparison, CP-EVP-FFT results are given in e, f for orientation #1

The equivalent Von Mises stress is clearly inhomogeneous in both phases as stress gradients develop in particles from the interface to the center of the particle, and in the matrix in regions between particles (hereinafter called channels) with the maximal value located at mid-distance between particles. The degree of heterogeneity appears stronger in the matrix. For both orientations, higher values of \( \sigma_{eq} \) are found in the regions between particles parallel to the tensile direction (Z-direction) and local stress hotspots are found close to sharp edges of the particles. For example, \( \sigma_{eq} \) is 1.5 times larger in the vicinity of particle corners than in the center of matrix channels for orientation #1 with the lowest \( l_{p} \). It is clearly observed that stress gradients are intensified in both phases with decreased interparticle spacing. Particle shape effects are stronger for lowest values of \( l_{p} \). In comparison, stress distribution obtained with CP-EVP-FFT shows that stress amplitudes are underestimated compared to MFDM-EVP-FFT. The analysis of statistical distribution of \( \sigma_{eq} \) reveals a smaller stress dispersion for CP-EVP-FFT, even if local hotspots are still found near particle edges, which is a common feature of CP-EVP-FFT and MFDM-EVP-FFT results.

Plots of the equivalent cumulated plastic strain \( \varepsilon_{eq}^{p} \) are reported in Fig. 5. Irrespective of interparticle spacing, it is observed that \( \varepsilon_{eq}^{p} \) is strongly inhomogeneous. Deformation bands develop in the ferrite matrix, and they are most pronounced for orientation #2 (Fig. 5c, d). For this case, it can be clearly seen that localization patterns conform to the spatial distribution of particles as bands develop along the longest uninterrupted paths into the microstructure. They are oriented at 45° from the tensile axis and the visible slip traces are consistent with the slip system with the highest Schmid factor in this single slip orientation. For orientation #1, bands also occur but plastic strain is more distributed around the particles, which is consistent with multiple activated slip systems on both possible plane families \( \left\{ {110} \right\} \) and \( \left\{ {211} \right\} \), see [64].

Fig. 5
figure 5

ad Spatial distributions of equivalent cumulated plastic strain \( \varepsilon_{eq}^{p} \) predicted with MFDM-EVP-FFT for a macroscopic strain \( E_{ZZ} = 0.2\% \), for two different particle interspacing distances \( l_{p} \) and two different crystallographic orientations of the matrix. 2-D sections containing the unit cell center are either parallel or perpendicular to the Z tensile direction (used axis systems are displayed in a). Each 2-D section represents the same physical length. For comparison, CP-EVP-FFT results are given in e, f for orientation #1 with a different scale range for better visibility

Plastic strain distribution undergoes significant change when \( l_{p} \) is decreased from \( l_{p} = 1.26\;\upmu{\text{m}} \) to \( l_{p} = 0. 10 5\;\upmu{\text{m}} \). It can be clearly seen that the region where very low levels of plastic strain are predicted grows from the matrix/particle interfaces to the mid-distance between particles. This results in a decrease of the magnitude of \( \varepsilon_{eq}^{p} \) for smaller interspacing distances. However, at the same time, localization bands become thinner and more intense as the region where plasticity occurs is further reduced in the channels between particles. This is confirmed by the analysis of the distribution of \( \varepsilon_{eq}^{p} /E_{eq}^{p} \), i.e. \( \varepsilon_{eq}^{p} \) normalized by the volume average of \( \varepsilon_{eq}^{p} \) over the matrix (not presented here), where longer tails are observed regarding high ratio values for low \( l_{p} \) values.

The spatial distribution of plastic strains predicted with CP-EVP-FFT is different from previous results, as it is more homogeneously distributed in the matrix, and plastic strain hotspots seem more tied to particle shape (Fig. 5e, f).

The spatial distribution of GND density \( \rho_{GND} \) obtained with MFDM-EVP-FFT is given in Fig. 6a–d. The 2-D sections show that \( \rho_{GND} \) increases from the center of matrix channels to matrix/particle interfaces and that the magnitude of \( \rho_{GND} \) is higher when particles are closer to each other. Distribution of GND densities is affected both by the shape of the particles and the crystallographic orientation of the matrix. For orientation #1, strong \( \rho_{GND} \) gradients developed along the interfaces as shown in both the parallel and perpendicular sections. Additionally, the parallel view shows that very high values of \( \rho_{GND} \) are found in the vicinity of particle corners and spread far into the matrix, demonstrating a strong shape effect on GND density for this orientation. This shape effect is not as pronounced in the case of orientation #2 (Fig. 6c, d). Layers of high GND density develop in the matrix close to particles but a network of low \( \rho_{GND} \) values is also formed (see arrows on Fig. 6c, d).

Fig. 6
figure 6

ad Spatial distributions of of GND densities \( \rho_{GND} \) predicted with MFDM-EVP-FFT for a macroscopic strain \( E_{ZZ} = 0.2\% \), for two different particle interspacing distances \( l_{p} \) and two different crystallographic orientations of the matrix. 2-D sections containing the unit cell center are either parallel or perpendicular to the Z tensile direction (used axis systems are displayed in a). Each 2-D section represents the same physical length. For comparison, CP-EVP-FFT results are given in e, f for orientation #1

For comparison, results obtained with CP-EVP-FFT given in Fig. 6e show that the post-calculated \( \rho_{GND} \) is localized at the particle’s interface/edges and no GND density pile-up develops and spreads inside the surrounding matrix. This is consistent with the earlier results reported in [45] for two-phase laminates. In addition, the magnitude of \( \rho_{GND} \) is 1–2 orders of magnitude lower than in the case of MFDM-EVP-FFT simulations.

Comparing the \( \sigma_{eq} \) and \( \rho_{GND} \) plots of Figs. 4 and 6, respectively, it can be seen that regions with stress hotspots correspond quite well to high GND densities, which indicates that significant plastic strain incompatibilities occur at these locations in the matrix. Moreover, the comparison between Figs. 5 and 6 shows that the spatial distribution of deformation bands at low strain is strongly affected by the consideration of non-local plasticity associated with the presence of GND pile-ups. Indeed, localization patterns, i.e. regions with high slip activity, conform exactly to regions with low GND. This is consistent with layers of GNDs obstructing the mean free path of mobile dislocations in the matrix. This mechanism cannot be described with conventional CP-EVP-FFT where GND densities are underestimated and not induced by pile-ups, and high plastic strain incompatibilities occur only at particle sharp edges. The effect of parameter \( \varvec{k}_{0} \) has not been studied in this paper but an increase of \( \varvec{k}_{0} \) could lead to higher local hardening and slightly higher GND densities close to phase interfaces as it was observed in previous works [45].

Interplay between particle interspacing and grain size effects

In this section, we investigate the interplay between particle interspacing and grain size effects using polycrystalline MMC microstructures. For that purpose, polycrystalline unit cells represented in Fig. 7 and discretized with 128 × 128 × 128 voxels were generated and used as input of MFDM-EVP-FFT numerical simulations. All unit cells are made of 8 cubic and randomly oriented grains (Fig. 7a shows the adopted morphology and orientations of the grains). A set of cubic grains was preferred over the use of Voronoi tessellation, to introduce a constant grain size instead of a distribution of grain sizes. Keeping the volume fraction of particles constant (close to 9.5%), regularly spaced hexagonal prismatic particles of different sizes were included in the polycrystalline unit cell. Same considerations as in “Effect of the particle interspacing distance” section were made about particle distribution and shape. Three unit cells were obtained by including 8, 64 or 216 particles (Fig. 7b). By considering two different voxel sizes for each of these three unit cells, two sets of results were obtained, respectively: (1) constant grain size of dg = 1.24 µm and interspacing distances of lp = 1 µm, 0.5 µm, and 0.33 µm, and (2) constant grain size of dg = 6.2 µm and interspacing distances of lp = 5 µm, 2.5 µm, and 1.66 µm.

Fig. 7
figure 7

Two-phase polycrystalline unit cells composed of eight cubic randomly oriented grains and 8, 64 or 216 regularly spaced particles and discretized with 128 × 128 × 128 voxels. Crystallographic orientations of grains are displayed in the standard stereographic triangle

Macroscopic responses

Macroscopic strain–stress curves of the two-phase polycrystalline unit cells are shown in Fig. 8. The change in grain size from dg = 6.2 µm to dg = 1.24 µm between the two sets of simulations clearly induces an increase of the macroscopic tensile stress. As it was reported in [46], this effect cannot be obtained with conventional CP-EVP-FFT, which is insensitive to grain size. Also, the dependence of macroscopic flow stress with grain size was shown to follow a Hall–Petch scaling law. However, no visible difference on the tensile curve is observed between the different \( l_{p} \) for a fixed grain size. It appears that the grain size effect is predominant on the macroscopic tensile response, at least for the investigated range of length scales. For comparison, results obtained with polycrystalline microstructures without addition of particles are reported for comparison as black dotted and plain lines, respectively for grain sizes dg = 6.2 µm and dg = 1.24 µm. The effective moduli of these microstructures being lower, a decrease of the slope in the elastic domain and a subsequent lower absolute macroscopic hardening are observed. However, the relative difference in macroscopic hardening between the two grain sizes is quite similar to that with particle-reinforced microstructures for the investigated strain range, confirming the apparent predominance of grain size effect.

Fig. 8
figure 8

Macroscopic stress–strain responses of the two-phase polycrystalline unit cells under uniaxial tension as predicted by MFDM-EVP-FFT, with grain size \( d_{g} = 1.24\;\upmu{\text{m}} \) and interspacing distances \( l_{p} \) ranging from 0.33 to 1 µm (solid lines), and \( d_{g} = 6.20\;\upmu{\text{m}} \) and interspacing distances \( l_{p} \) ranging from 1.55 to 5 µm (dotted lines). Black dotted and plain lines give the tensile response of polycrystalline unit cells without particles, with respective grain sizes of \( {d}_{{g}} = 6.20\;\upmu{\text{m}} \) and \( {d}_{{g}} = 1.24\;\upmu{\text{m}} \)

Analysis of intra-granular mechanical fields

The scalar micromechanical fields defined in “Outputs” section, predicted at \( E_{ZZ} = 0. 2\% \), are analyzed here. Only results for dg = 1.24 µm are presented. Similar trends were obtained for dg = 6.2 µm.

The field of equivalent Von Mises stress \( \sigma_{eq} \) inside four adjacent grains is shown in Fig. 9a, b. As was observed in the cases of single crystal matrix in “Effect of the particle interspacing distance” section, stress is strongly inhomogeneous. Stress heterogeneity is observed comparing the fields inside the ferrite phase and elastic particles, and now also comparing the fields inside different grains. Two kinds of stress hotspots are identified in the grains:

Fig. 9
figure 9

Spatial distribution of equivalent Von Mises stress \( \sigma_{eq} \) for \( E_{ZZ} = 0.2\% \) predicted with MFDM-EVP-FFT for unit cells with \( d_{g} = 1.24\;\upmu{\text{m}} \) and a 8 or b 216 particles. 2-D sections through the center of the unit cell and in the middle of grains #1 to #4

  1. i.

    The highest stress values are found at grain boundaries near to a particle (see white boxes in Fig. 9b). These hotspots are present all along grain boundaries and are separated from each other by regions of low \( \sigma_{eq} \).

  2. ii.

    Local maxima of \( \sigma_{eq} \) are also observed between particles in the Z-direction, i.e. the tensile direction (see yellow boxes in Fig. 9b).

With the decrease of particle spacing, locations described in (i) become frequent along the grain boundaries with a slight reduction of the highest value of \( \sigma_{eq} \).

The fields of the equivalent plastic strain \( \varepsilon_{eq}^{p} \) are reported in Fig. 10a, b. For the single slip orientations among the 8 considered grains, deformation bands oriented at 45° from principal axis of tensile loading are clearly visible and can be linked to best oriented slip systems (as it was observed for orientation #2 in “Effect of interspacing distance on intra-granular mechanical fields” section). Interestingly, plastic strain values are higher in matrix channels in the center of grains, away from grain boundaries, which act as new impenetrable obstacles to slip. It is important to note that because of the cubic distribution of particles chosen here (no clustering), this statistical homogeneous decrease of particle interspacing distance results in the multiplication of possible channel paths for mobile dislocations, as it is observed from the increased number of slip lines between lp = 1 µm and lp = 0.33 µm.

Fig. 10
figure 10

Spatial distribution of equivalent plastic strain \( \varepsilon_{eq}^{p} \) for \( E_{ZZ} = 0.2\% \) predicted with MFDM-EVP-FFT for unit cells with \( d_{g} = 1.24\;\upmu{\text{m}} \) and a 8 or b 216 particles. 2-D sections through the center of the unit cell and in the middle of grains #1 to #4

Finally, plots of \( \rho_{GND} \) are reported on Fig. 11a, b. GND pile-ups develop at both matrix/particle interfaces and at grain boundaries. However, they are larger and more intense close to grain boundaries, especially in grain #3. Furthermore, layers of highest GND density are not distributed along grain boundaries. Indeed, regions where higher values of \( \rho_{GND} \) are obtained correspond to the stress hotspots previously described, indicating local high plastic strain incompatibilities. Interestingly, stress and GND density hotspots are located at the termination of plastic strain bands observed in Fig. 10b.

Fig. 11
figure 11

a, b Spatial distribution of GND density \( \rho_{GND} \) for \( E_{ZZ} = 0.2\% \) predicted with MFDM-EVP-FFT for unit cells with \( d_{g} = 1.24\;\upmu{\text{m}} \) and a 8 or b 216 particles. 2-D sections through the center of the unit cell and in the middle of grains #1 to #4. c Profiles of intra-granular GND density \( \rho_{GND} \) along different paths indicated in a, b by colored arrows

When particles become closer to each other for smaller \( l_{p} \), the \( \rho_{GND} \) concentrations progressively disappear in grains #1, #2 and #4 of Fig. 11b. The careful analysis of profiles of \( \rho_{GND} \) along specific paths in grain #3, illustrated in Fig. 11c, also indicates a slight decrease of maximum values of \( \rho_{GND} \) at grain boundaries. This can be explained by the following mechanism: for lp = 0.33 µm, the number of observed deformation bands carrying mobile dislocations is higher, but more particles are intersected by dislocations on each of these paths as the interspacing distance is smaller. Thus, as overall deformation increases, more intersections occur, where dislocations can be trapped in the vicinity of particles, increasing the local \( \rho_{GND} \) magnitude around particles. Therefore, a lower amount of dislocations contributes to the intensification of GND pile-ups at grain boundaries. Concomitantly, a more GND pile-ups can be formed at grain boundaries as more intersections with deformation paths are observed. These two phenomena are likely to balance each other out for the ranges of length scales investigated here, justifying that no interparticle spacing effect was observed on the macroscopic tensile response in “Macroscopic responses” section for a given grain size.

Considering that, as particles get even closer to each other, and below a critical particle interspacing distance to grain size ratio, the GND pile-ups at matrix/particle interfaces would be more detrimental to plastic strain than pile-ups at grain boundaries, then the particle interspacing effect should take over the grain size effect. This, however, is not observed here, as particle interspacing distances are probably still too high relatively to grain size. More importantly, observations made here are specific, and maybe restricted, to the chosen uniform distribution of particles with cubic periodicity. For clustered distributions of particles inside grains, local distances between neighboring particles could be very low compared to grain size even for a large mean particle interspacing. The presence of these particle clusters are likely to be detrimental to the emergence of deformation patterns, and could lead to early and rapid damage initiation in the vicinity of particles during loading, as was recently shown in [65]. MFDM-EVP-FFT simulations on synthetic microstructures with particle clusters, to study the effect of non-uniform particle distributions on the mechanical behavior of Fe–TiB2 will be presented in a future contribution.

Conclusions

A spectral formulation called MFDM-EVP-FFT that includes GND and SSD effects through Field Dislocation Mechanics as an extension of the EVP-FFT formulation [39, 45] was used in this paper. This FFT-based approach was able to successfully describe internal length scales effects in a Fe–TiB2 metal matrix composite considering both elastic particles and elasto-viscoplastic matrix. The matrix was considered to be either single crystalline (periodic unit cell with one ferrite grain without grain boundaries) or polycrystalline (periodic unit cell with eight grains including the presence of grain boundaries). Numerical simulations with different voxelized synthetic unit cells obtained with a periodic cubic distribution of particles and different particle interspacing distances were performed with 64 × 64 × 64 voxels or 128 × 128 × 128 voxels.

Numerical results for single-crystalline matrix demonstrated a strong dependence of equivalent stress, cumulated plastic strain and GND density with the particle interspacing distance. In the non-local formulation of MFDM-EVP-FFT, dependence of mechanical behavior and distribution of intra-granular mechanical fields is related to the formation of higher GND density pile-ups at matrix/particle interfaces as opposed to conventional crystal plasticity (CP-EVP-FFT). MFDM-EVP-FFT results show spatial correlations between low and high values of GND density and respectively high and low values of equivalent cumulated plastic strain. Indeed, numerical simulations were able to describe that, with decreased distance between neighboring particles, GND pile-ups spread further into grain interiors. Thus, regions in matrix channels where slip activity is inhibited become larger leading to a more localized network of deformation bands and stronger slip gradients inside grains. Once GND pile-ups build up around particles, slip is further constrained in the matrix due to these GND pile-ups that conform with particle location and shape. This mechanism induces stronger hardening in a similar way to mechanism described for example in [66,67,68]. This effect was shown to depend strongly on the shape of particles and the crystallographic orientation in the matrix. Indeed, for orientation #2 considered in “Effect of the particle interspacing distance” section, it was seen that GND accumulation close to particles does not significantly impede the single slip activity as regions of low \( \rho_{GND} \) conform to the activated slip system and the longest available paths for dislocation motion in the microstructure, irrespective of interparticle spacing. Consequently, the hardening effect on the macroscopic behavior induced by the formation of GND pile-ups is less pronounced for orientation #2 than for orientation #1, as reported on Fig. 3.

Following the analysis of numerical results where a polycrystalline matrix is described, grain-size effect was shown to be predominant over particle interspacing effect for the range of investigated length scales and it was attributed to a too low grain size to particle spacing ratio and the choice of a periodic cubic distribution of particles in grain interiors. Higher GND densities and more intense deformation patterns are likely to be observed in clustered configurations with non-uniform particle distribution. The relationship between the degree of clustering, the macroscopic hardening and the spatial distribution will be discussed in a future contribution.

Availability of data and materials

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

References

  1. Hall EO. The deformation of ageing and mild steels. Proc Phys Soc Lond B. 1951;64:747–53.

    Article  Google Scholar 

  2. Petch NJ. The cleavage strength of polycrystals. J Iron Steel Inst. 1953;174:25–8.

    Google Scholar 

  3. Ashby MF. Deformation of plastically non-homogeneous materials. Philos Mag. 1970;21:399–424.

    Article  Google Scholar 

  4. Cherkaoui M, Capolungo L. Atomistic and continuum modeling of nanocrystalline materials: deformation mechanisms and scale transition. Berlin: Springer; 2009.

    Book  Google Scholar 

  5. Lloyd DJ. Particle reinforced aluminum and magnesium matrix composites. Int Mater Rev. 1994;39:1–23.

    Article  Google Scholar 

  6. Nan CW, Clarke DR. The influence of particle size and particle fracture on the elastic/plastic deformation of metal matrix composites. Acta Mater. 1996;44(9):3801–11.

    Article  Google Scholar 

  7. Fleck NA, Muller GM, Ashby MF, Hutchinson JW. Strain gradient plasticity: theory and experiment. Acta Metall Mater. 1994;42:475–87.

    Article  Google Scholar 

  8. Fleck NA, Hutchinson JW. Strain gradient plasticity. In: Hutchinson JW, Wu TY, editors. Advances in applied mechanics, vol. 33. New York: Academic Press; 1997. p. 295–361.

    Google Scholar 

  9. Fleck NA, Hutchinson JW. A reformulation of strain gradient plasticity. J Mech Phys Solids. 2001;49(10):2245–71.

    Article  MATH  Google Scholar 

  10. Aifantis EC. On the microstructural origin of certain inelastic models. J Eng Mater Technol. 1984;106:326–30.

    Article  Google Scholar 

  11. Aifantis EC. The physics of plastic deformation. Int J Plast. 1987;3:211–47.

    Article  MATH  Google Scholar 

  12. Zbib HM, Aifantis EC. On the localization and postlocalization behavior of plastic deformation—part I: on the initiation of shear bands. Res Mech. 1988;23:261–77.

    Google Scholar 

  13. Zbib HM, Aifantis EC. A gradient-dependent theory of plasticity: application to metal and soil instabilities. Appl Mech Rev. 1989;42:295–304.

    Article  MATH  Google Scholar 

  14. Acharya A, Bassani JL. Lattice incompatibility and a gradient theory of crystal plasticity. J Mech Phys Solids. 2000;48(8):1565–95.

    Article  MathSciNet  MATH  Google Scholar 

  15. Gao H, Huang Y, Nix WD, Hutchinson JW. Mechanism-based strain gradient plasticity-I. J Mech Phys Solids. 1999;47:1239–63 [see also Huang Y, Gao H, Nix WD, Hutchinson JW. Mechanism-based strain gradient plasticity-II. J Mech Phys Solids. 2000;48:99–128].

  16. Fleck NA, Hutchinson JW. A phenomenological theory for strain gradient effects in plasticity. J Mech Phys Solids. 1993;41(12):1825–57.

    Article  MathSciNet  MATH  Google Scholar 

  17. Acharya A, Beaudoin AJ. Grain-size effect in viscoplastic polycrystals at moderate strains. J Mech Phys Solids. 2000;48(10):2213–30.

    Article  MATH  Google Scholar 

  18. Wulfinghoff S, Böhlke T. Gradient crystal plasticity including dislocation-based work-hardening and dislocation transport. Int J Plast. 2015;69:152–69.

    Article  Google Scholar 

  19. Kröner E. Continuum theory of defects. In: Balian R, et al., editors. Physics of defects. Les Houches Session 35 North Holland: New York; 1981. p. 215–315.

    Google Scholar 

  20. Acharya A. A model of crystal plasticity based on the theory of continuously distributed dislocations. J Mech Phys Solids. 2001;49:761–84.

    Article  MATH  Google Scholar 

  21. Acharya A. Driving forces and boundary conditions in continuum dislocation mechanics. Proc R Soc Lond A. 2003;459:1343–63.

    Article  MathSciNet  MATH  Google Scholar 

  22. Acharya A. Constitutive analysis of finite deformation field dislocation mechanics. J Mech Phys Solids. 2004;52:301–16.

    Article  MathSciNet  MATH  Google Scholar 

  23. Nye JF. Some geometrical relations in dislocated crystals. Acta Mater. 1953;1:153–62.

    Article  Google Scholar 

  24. Acharya A, Roy A. Size effects and idealized dislocation microstructure at small scales: predictions of a phenomenological model of mesoscopic field dislocation mechanics: part I. J Mech Phys Solids. 2006;54:1687–710.

    Article  MathSciNet  MATH  Google Scholar 

  25. Roy A, Acharya A. Size effects and idealized dislocation microstructure at small scales: predictions of a Phenomenological model of mesoscopic field dislocation mechanics: part II. J Mech Phys Solids. 2006;54:1711–43.

    Article  MathSciNet  MATH  Google Scholar 

  26. Roy A, Acharya A. Finite element approximation of field dislocation mechanics. J Mech Phys Solids. 2005;53:143–70.

    Article  MATH  Google Scholar 

  27. Taupin V, Varadhan S, Chevy J, Fressengeas C, Beaudoin AJ, Montagnat M, Duval P. Effects of size on the dynamics of dislocations in ice single crystals. Phys Rev Lett. 2007;99:155507.

    Article  Google Scholar 

  28. Puri S, Roy A, Acharya A, Dimiduk D. Modeling dislocations sources and size effects at initial yield in continuum plasticity. J Mech Mater Struct. 2009;4(9):1603–18.

    Article  Google Scholar 

  29. Puri S, Das A, Acharya A. Mechanical response of multicrystalline thin films in mesoscale field dislocation mechanics. J Mech Phys Solids. 2011;59:2400–17.

    Article  MATH  Google Scholar 

  30. Puri S, Roy A. Plastic deformation of multicrystalline thin films: grain size distribution vs. grain orientation. Comput Mater Sci. 2012;52:20–4.

    Article  Google Scholar 

  31. Taupin V, Varadhan S, Fressengeas C, Beaudoin AJ. Directionality of yield point in strain-aged steels: the role of polar dislocations. Acta Mater. 2008;56:3002–10.

    Article  Google Scholar 

  32. Roy A, Puri S, Acharya A. Phenomenological mesoscopic field dislocation mechanics, lower-order gradient plasticity and transport of mean excess dislocation density. Model Simul Mater Sci Eng. 2007;15:167–80.

    Article  Google Scholar 

  33. Richeton T, Wang GF, Fressengeas C. Continuity constraints at interfaces and their consequences on the work hardening of metal-matrix composites. J Mech Phys Solids. 2011;59:2023–43.

    Article  MathSciNet  MATH  Google Scholar 

  34. Moulinec H, Suquet P. A fast numerical method for computing the linear and non linear properties of composites. C R de l’Académie des Sci Paris II. 1994;318:1417–23.

    MATH  Google Scholar 

  35. Moulinec H, Suquet P. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput Methods Appl Mech Eng. 1998;157:69–94.

    Article  MathSciNet  MATH  Google Scholar 

  36. Eyre DJ, Milton GW. A fast numerical scheme for computing the response of composite using grid refinement. Eur Phys J Appl Phys. 1999;6:41–7.

    Article  Google Scholar 

  37. Michel JC, Moulinec H, Suquet P. A computational scheme for linear and non-linear composites with arbitrary phase contrast. Int J Numer Methods Eng. 2001;52:139–60.

    Article  Google Scholar 

  38. Lebensohn RA. N-site modeling of a 3D viscoplastic polycrystal using fast-Fourier transform. Acta Mater. 2001;49:2723–37.

    Article  Google Scholar 

  39. Lebensohn RA, Kanjarla AK, Eisenlohr P. An elasto-viscoplastic formulation based on Fast Fourier Transforms for the prediction of micromechanical fields in polycrystalline materials. Int J Plast. 2012;32–33:59–69.

    Article  Google Scholar 

  40. Suquet P, Moulinec H, Castelnau O, Montagnat M, Lahellec N, Grennerat F, Duval P, Brenner R. Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep. Procedia IUTAM. 2012;3:76–90.

    Article  Google Scholar 

  41. Brenner R, Beaudoin AJ, Suquet P, Acharya A. Numerical implementation of static Field Dislocation Mechanics theory for periodic media. Philos Mag. 2014;94(16):1–24.

    Article  Google Scholar 

  42. Berbenni S, Taupin V, Djaka KS, Fressengeas C. A numerical spectral approach for solving elasto-static field dislocation and g-disclination mechanics. Int J Solids Struct. 2014;51:4157–75.

    Article  Google Scholar 

  43. Djaka KS, Villani A, Taupin V, Capolungo L, Berbenni S. Field dislocation mechanics for heterogeneous elastic materials: a numerical spectral approach. Comput Methods Appl Mech Eng. 2017;315:921–42.

    Article  MathSciNet  Google Scholar 

  44. Djaka KS, Taupin V, Berbenni S, Fressengeas C. A numerical spectral approach to solve the dislocation density transport equation. Model Simul Mater Sci Eng. 2015;23:065008(27 pp).

    Article  Google Scholar 

  45. Djaka KS, Berbenni S, Taupin V, Lebensohn RA. A FFT-based numerical implementation of mesoscale field dislocation mechanics: application to two-phase laminates. Int J Solids Struct. 2019. https://doi.org/10.1016/j.ijsolstr.2018.12.027.

    Article  Google Scholar 

  46. Berbenni S, Taupin V, Lebensohn RA. A fast-Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals. J Mech Phys Solids. 2020;135:103808.

    Article  MathSciNet  Google Scholar 

  47. Belkhabbaz A, Brenner R, Rupin N, Bacroix B, Fonseca J. Prediction of the overall behavior of a 3D microstructure of austenitic steel by using FFT numerical scheme. Proc Eng. 2011;10:1883–8.

    Article  Google Scholar 

  48. Rovinelli A, Proudhon H, Lebensohn RA, Sangid MD. Assessing the reliability of fast Fourier transform-based crystal plasticity simulations of a polycrystalline material near a crack tip. Int J Sol Struct. 2020;184:153–66.

    Article  Google Scholar 

  49. ARCELOR Research group. Patent EP 1 897 963 A1, Bulletin 2008/11, 20p.

  50. Lartigue-Korinek S, Walls M, Haneche N, Cha L, Mazerolles L, Bonnet F. Interfaces and defects in successfully hot-rolled steel-based composite Fe–TiB2. Acta Mater. 2015;98:297–305.

    Article  Google Scholar 

  51. Hadjem-Hamouche Z, Chevalier JP, Cui Y, Bonnet F. Deformation behavior and damage evolution in a new titanium diboride (TiB2) steel-based composite. Steel Res Int. 2012;83:538–45.

    Article  Google Scholar 

  52. Dammak M, Gaspérini M, Barbier D. Microstructural evolution of iron-based metal-matrix composites submitted to simple shear. Mater Sci Eng A. 2014;616:123–31.

    Article  Google Scholar 

  53. Gaspérini M, Dammak M, Franciosi P. Stress estimates for particle damage in Fe–TiB2 metal matrix composites from experimental data and simulation. Eur J Mech A/Solids. 2017;64:85–98.

    Article  Google Scholar 

  54. Mura T. Continuous distribution of moving dislocations. Philos Mag. 1963;89:843–57.

    Article  Google Scholar 

  55. Mura T. Periodic distribution of dislocations. Proc R Soc Lond A. 1964;280:528–44.

    Article  MathSciNet  MATH  Google Scholar 

  56. Mecking H, Kocks UF. Kinetics of flow and strain-hardening. Acta Metall. 1981;29:1865–75.

    Article  Google Scholar 

  57. Kocks UF. Laws for work-hardening and low-temperature creep. ASME J Eng Mater Technol. 1976;98:76–85.

    Article  Google Scholar 

  58. Mura T. Micromechanics of defects in solids. Dordrecht: Edition Kluwer Academic Publishers; 1987.

    Book  MATH  Google Scholar 

  59. Willot F. Fourier-based schemes for computing the mechanical response of composites with accurate local fields. C R Mécanique. 2015;343:232–45.

    Article  Google Scholar 

  60. Okamoto NL, Kusakari M, Tanaka K, Inui H, Otani S. Anisotropic elastic constants and thermal expansivities in monocrystal CrB2, TiB2 and ZrB2. Acta Mater. 2010;58(1):76–84.

    Article  Google Scholar 

  61. Groeber MA, Jackson MA. DREAM.3D: a digital representation environment for the analysis of microstructure in 3D. Integr Mater. 2014;3(1):56–72.

    Google Scholar 

  62. Tomita Y, Higa Y, Fujimoto T. Modeling and estimation of deformation behavior of particle-reinforced metal-matrix composite. Int J Mech Sci. 2000;42(11):2249–60.

    Article  MATH  Google Scholar 

  63. Zhu HT, Zbib HM, Aifantis EC. Strain gradients and continuum modeling of size effect in metal matrix composites. Acta Mech. 1997;121(1–4):165–76.

    Article  MATH  Google Scholar 

  64. Weinberger CR, Loyce BL, Battaile CC. Slip planes in bcc transition metals. Int Mater Rev. 2013;58:296–314.

    Article  Google Scholar 

  65. Wang D, Shanthraj P, Springer H, Raabe D. Particle-induced damage in Fe–TiB2 high stiffness metal matrix composites steels. Mater Des. 2018;160:557–71.

    Article  Google Scholar 

  66. Cleveringa HHM, Van Der Giessen E, Needleman Z. Comparison of discrete dislocation and continuum plasticity predictions for a composite material. Acta Mater. 1997;45(8):3163–79.

    Article  Google Scholar 

  67. Bassani JL, Needleman A, Van Der Giessen E. Plastic flow in a composite: a comparison of nonlocal continuum and discrete dislocation predictions. Int J Solids Struct. 2001;38:833–53.

    Article  MATH  Google Scholar 

  68. Fredriksson P, Gudmundson P, Mikkelsen LP. Finite element implementation and numerical issues of strain gradient plasticity with application to metal matrix composites. Int J Solids Struct. 2009;46:3977–87.

    Article  MATH  Google Scholar 

Download references

Funding

This study was conducted in the frame of a Carnot ARTS project in collaboration with J.-P. Chevalier and co-workers at PIMM Laboratory (CNAM, Arts et Métiers Paris Tech, CNRS). The authors thank the French State (ANR) through the program “Investment in Future” (LabEx “DAMAS” referenced as ANR-11-LABX-0008-01), ArcelorMittal Research Development and the Grand Est region for financial support. RAL acknowledges support from Los Alamos National Laboratory’s Laboratory-Directed Research and Development (LDRD) program.

Author information

Authors and Affiliations

Authors

Contributions

The MFDM-EVP-FFT model was initially developed by SB as an extension of the EVP-FFT formulated by RAL. JG conceived the application of MFDM-EVP-FFT to the studied MMC material, collected and interpreted the numerical results, and was the major contributor in writing the manuscript. SB, NG and RAL provided valuable suggestions. All authors read and approved the final manuscript.

Corresponding author

Correspondence to J. Genée.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Genée, J., Berbenni, S., Gey, N. et al. Particle interspacing effects on the mechanical behavior of a Fe–TiB2 metal matrix composite using FFT-based mesoscopic field dislocation mechanics. Adv. Model. and Simul. in Eng. Sci. 7, 6 (2020). https://doi.org/10.1186/s40323-020-0141-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-020-0141-z

Keywords