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

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.

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 singlecrystal 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 spacetime 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 TiB 2 (hexagonal compact (HCP) crystallographic structure). This new composite is attractive for industrial applications due to the very high elastic moduli and low density of TiB 2 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-TiB 2 [51,52]. It was demonstrated that cracking of TiB 2 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-TiB 2 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 elastoviscoplastic 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 TiB 2 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.

MFDM equations
To solve the displacement u and stress σ 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 ∪ S u ): where 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 U p , which results from dislocation motion, and the average elastic (or lattice) distortion U 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 α . The mesoscale FDM theory is based on an averaged value of the tensor α. 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: The mobility of SSDs is represented by the mesoscale plastic distortion rate L p , where the averaging procedure was defined in [24]. The space-time evolution of the average dislocation density tensor α is obtained from the conservation of dislocation flux and is prescribed as [54,55]: Constitutive specifications on the dislocation velocity v and on the slip distortion rate L p are given from thermodynamic considerations following the theory introduced by Acharya and Roy [24]. Furthermore, plastic flow incompressibility is considered. The GND velocity v is prescribed as follows: where g is the glide force parallel to v and v is the magnitude of 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: where N is the total number of slip systems ( N = 24 for the BCC structure in ferrite: α-Fe), γ s is the slip rate on slip system s , ζ is a material constant, b is the magnitude of the Burgers vector, τ c is the shear strength and µ is the isotropic elastic shear modulus of the material.
In crystal plasticity, the plastic distortion rate tensor L p due to slip is defined as: where b s and n s denote, respectively for each slip system s , the slip direction and the slip plane unit normal. The constitutive equation for γ s introduced in Eqs. 5 and 6 is given by a classic viscoplastic flow rule as a power law: where m is the strain rate sensitivity of the material, τ s is the resolved shear stress, γ 0 is the reference slip rate and τ 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: The evolution law for the shear strength τ c follows the same hypotheses as the strainhardening model developed by [29], which is an extension of earlier models derived by [17,56]: where τ 0 is the yield strength due to lattice friction, τ s is the saturation stress, θ 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. α = 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: where ξ is the Fourier vector, Ŵ 0 is the modified Green tensor associated with the homogeneous elastic moduli C 0 , and τ ij = σ ij − C 0 ijkl ε kl represents the stress polarization tensor field due to heterogeneities. As for notation, ε(ξ ) is the continuous Fourier transform of tensor field ε(x) . 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 (n) are approximated by auxiliary fields (n) ij and e (n) kl respectively. The stress polarization tensor then becomes: The augmented Lagrangian scheme also requires the nullification of the residual R , which depends on the stress and strain tensors σ (n+1) and ε (n+1) (Eq. 12). This non-linear equation is solved using a Newton-Raphson iterative scheme.
Once the convergence is reached for σ (n+1) and ε (n+1) , the new guess for the auxiliary stress field is given by using Uzawa's descent algorithm: The iterative algorithm is stopped when the normalized average differences between the stress fields σ and and the strain fields ε and 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 E = ε (n) is applied to the periodic unit cell V in the form of: In cases of mixed boundary conditions with imposed macroscopic strain rate Ė ij and stress Σ ij , the (n + 1)-guess of the macroscopic strain E (n+1) ij 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 α t+�t ij where 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 η as: where β = −log(ε M ) and 2p correspond to the high-frequency damping parameter and the order of the filter respectively. Applying the filter to Eq. 3 leads to: 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 t is fixed according to a Courant-Friedrichs-Lewy (CFL) condition: where t CFL is the refined time step, δ is the voxel size, c is a user-specified constant and v max is the maximal GND velocity.

Description of material
The material studied here is a ferritic-steel-based composite reinforced with a 9.5% volume fraction of TiB 2 particles. The composite was produced by ArcelorMittal Research SA, by ingot casting with in situ precipitation of the ceramic TiB 2 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 TiB 2 particles are uniformly distributed in the matrix, elongated preferentially in the rolling direction and exhibit hexagonal-like sections. Moreover, the crystal structure of the TiB 2 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 TiB 2 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-TiB 2 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 TiB 2 , respectively. For the elasto-viscoplatic matrix, the material parameters related to slip rule and GND velocity ( γ 0 , m in Eq. 7 and ζ 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 ( τ 0 , τ s and θ 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 × 10 −10 m.
Given the body-centered cubic structure of the α-Fe lattice, both {110} 111 and {211} 111 slip systems are considered (i.e. a total of 24 slip systems, with {123} 111 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 Ė ZZ = 0.001 s −1 . No specific condition was considered   Table 2 Material parameters related to plasticity used for numerical simulationṡ 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, ε 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 ( τ 0 , τ s and θ 0 described in "MFDM equations" section). To this end, a periodic unit cell representative of the Fe-TiB 2 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 d g = 4.5 µm and standard deviation of about 0.4 × d g . TiB 2 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 d p = 1.7 µm and standard deviation of about 0.6 × d g . Finally, the c-axis of the HCP lattice is oriented along Z-axis for all particles.
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%: i. The equivalent Von Mises stress σ eq , in order to evidence stress hotspots in the microstructures: where ε p is the symmetric part of U p . iii. The L 2 -norm α of the Nye tensor or alternatively the GND density, defined as: The simplified analysis of α 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 = s y = 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: 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.105 µm to l p = 1.26 µ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 111 lies parallel to the tensile direction Z. For the second, single slip is the preferred slip mode (with 112 111 system exhibiting the highest Schmid factor) and 5, 6, 27 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 µm down to l p = 0.105 µ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 (21)  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]. 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 σ 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 σ 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. 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 σ 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, σ 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 σ 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 ε p eq are reported in Fig. 5. Irrespective of interparticle spacing, it is observed that ε p eq 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 {110} and {211} , see [64].
Plastic strain distribution undergoes significant change when l p is decreased from l p = 1.26 µm to l p = 0.105 µ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 middistance between particles. This results in a decrease of the magnitude of ε p eq 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 ε p eq /E p eq , i.e. ε p eq normalized by the volume average of ε p eq 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 ρ GND obtained with MFDM-EVP-FFT is given in Fig. 6a-d. The 2-D sections show that ρ GND increases from the center of matrix channels to matrix/particle interfaces and that the magnitude of ρ 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 ρ 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 ρ 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 ρ GND values is also formed (see arrows on Fig. 6c, d).
For comparison, results obtained with CP-EVP-FFT given in Fig. 6e show that the post-calculated ρ 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 ρ GND is 1-2 orders of magnitude lower than in the case of MFDM-EVP-FFT simulations. Comparing the σ eq and ρ 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 k 0 has not been studied in this paper but an increase of 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 d g = 1.24 µm and interspacing distances of l p = 1 µm, 0.5 µm, and 0.33 µm, and (2) constant grain size of d g = 6.2 µm and interspacing distances of l p = 5 µm, 2.5 µm, and 1.66 µm.

Macroscopic responses
Macroscopic strain-stress curves of the two-phase polycrystalline unit cells are shown in Fig. 8. The change in grain size from d g = 6.2 µm to d g = 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 d g = 6.2 µm and d g = 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 Fig. 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 similar to that with particle-reinforced microstructures for the investigated strain range, confirming the apparent predominance of grain size effect.

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 d g = 1.24 µm are presented. Similar trends were obtained for d g = 6.2 µm.
The field of equivalent Von Mises stress σ 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: 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 σ eq . ii. Local maxima of σ 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 σ eq . The fields of the equivalent plastic strain ε p eq 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 l p = 1 µm and l p = 0.33 µm.
Finally, plots of ρ 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 ρ 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.
When particles become closer to each other for smaller l p , the ρ GND concentrations progressively disappear in grains #1, #2 and #4 of Fig. 11b. The careful analysis of profiles of ρ GND along specific paths in grain #3, illustrated in Fig. 11c, also indicates a slight decrease of maximum values of ρ GND at grain boundaries. This can be explained by the following mechanism: for l p = 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 ρ 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-TiB 2 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-TiB 2 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 Fig. 11 a, b Spatial distribution of GND density ρ GND for E ZZ = 0.2% predicted with MFDM-EVP-FFT for unit cells with d g = 1.24 µ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 ρ GND along different paths indicated in a, b by colored arrows 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 ρ 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.