The effects of heterogeneous mechanical properties on the response of a ductile material

We investigate numerically the small-strain, elastic–plastic response of statistically isotropic materials with non-uniform spatial distributions of mechanical properties. The numerical predictions are compared to simple bounds derived analytically. We explore systematically the effects of heterogeneity on the macroscopic stiffness, strength, asymmetry, stability and size dependence. Monte Carlo analyses of the response of statistical volume elements are conducted at different strain triaxiality using computational homogenisation, and allow exploring the macroscopic yield behaviour of the heterogeneous material. We illustrate quantitatively how the pressure-sensitivity of the yield surface of the solid increases with heterogeneity in the elastic response. We use the simple analytical models developed here to derive an approximate scaling law linking the fatigue endurance threshold of metallic alloys to their stiffness, yield strength and tensile strength.


Analytical predictions of the uniaxial response
We begin by deriving approximate analytical bounds to the macroscopic uniaxial stress-strain response of elastic-plastic solids with microscopic variations of either the Young's modulus E, the yield strength σ y or the hardening modulus H, defined as the slope of the true stress-strain curve in the inelastic phase; the elastic Poisson's ratio ν is taken as homogeneous and equal to 0.3. We stress here that in real materials the degrees of heterogeneity in these three local properties might be correlated to some extent; here however we do not aim at representing a particular material, but rather we focus on the effects of local variations of each of these three properties separately.
We assume that the material consists of a cuboidal array of N elastic-plastic homogeneous cubic SVEs (or 'cells') with linear strain hardening; the relevant microscopic material properties are assumed to follow a uniform random distribution as follows with i = 1 . . . N . Here E max , E min , σ max y , σ min y , H max and H min are maximum and minimum values associated with the ranges of variation of the three mechanical properties; r i denotes uniformly distributed, uncorrelated random numbers in the interval [0,1]; the average values of the mechanical properties are denoted as E 0 , σ y0 and H 0 . The choice of an uncorrelated, uniform probability density for the mechanical properties of each phase of the composite was based on the ease of implementation, on the simplicity of its analytical treatment, and the fact that it ensures the maximum possible variance for properties defined in a finite interval, such to enhance the effects of heterogeneity investigated here; we will show below that predictions are approximately insensitive to the shape of the (symmetric) probability density function considered but only to its variance.
The expected value of an upper bound for the macroscopic stress-strain response of the heterogeneous solid can be obtained assuming equal strain in all cells (Voigt bound) while an expected value of a lower bound can be calculated assuming equal stress in all cells (Reuss bound). In the following we will assume for simplicity that all cells are in a state of uniaxial stress during deformation.
Heterogeneity in Young's modulus. First we consider the case where the Young's modulus E (stiffness of the SVEs) varies randomly across the solid according to Eq. (1), while the yield strength and hardening modulus remain uniform, i.e., �σ y = 0 and H = 0 . We also assume H 0 = 0 , i.e., all cells are perfectly plastic.
Expected value of upper bound. The stress-strain curves of each cell of the statistical mesh are sketched in Fig. 1a; it is assumed that all cells undergo the same uniaxial microscopic strain, equal to the current macroscopic strain ε . For ε < ε min , where ε min = σ y0 /E max , all cells respond elastically; if ε min ≤ ε ≤ ε max , where ε max = σ y0 /E min , some of the cells respond elastically while others undergo plasticity; all cells deform plastically if ε > ε max .
For a given macroscopic applied strain ε , a number N p of cells is in the plastic regime; we define the fraction of plastic cells as (1) E i (r i ) = E min + r i (E max − E min ) = E min + r i �E where N e is the number of elastically deforming cells. If ε < ε min , all cells undergo a purely elastic response, P p = 0 . If ε > ε max , P p = 1 . In the transition region, ε min ≤ ε ≤ ε max , with reference to Fig. 1a, the fraction of plastic cells is given by where the variables a and b represent the length of the two segments sketched in the figure. In this phase, plastic and elastic regions coexist. The homogenized (or macroscopic) stress σ is obtained by averaging the stresses in all cells, as follows Combining Eqs. (4)- (7), and recalling that r i are uniformly distributed random variables, the expected value of the homogenized stress above can be calculated as  For ε > ε max , P p = 1 and σ = σ y0 .
Expected value of lower bound. Now assume that all cells are connected in series and experience an increasing but uniform macroscopic stress σ . If σ < σ y0 , the homogenized engineering strain ε is given by assuming small deformations. Recalling that r i are uniformly distributed random variables, we can obtain the expected value of the lower bound as which defines the macroscopic elastic behaviour. The response is linear elastic for σ < σ y0 ; when the applied stress is equal to the yield stress, the strain is undetermined (for this perfectly plastic material) and σ = σ y0 .
Heterogeneity in yield stress. Next, we consider the case of the local yield stress σ y varying randomly across the material domain according to Eq. (2), while Young's modulus and hardening modulus remain uniform, i.e. E = 0 , H 0 = 0 and H = 0.
Expected value of upper bound. We now assume that all cells are subject to the uniform macroscopic applied uniaxial strain ε . With reference to Fig. 1b, three regions of the response can be identified; for ε < ε min , where ε min = σ min y /E 0 , all cells respond elastically; a transition region exists for ε min ≤ ε ≤ ε max = σ max y /E 0 , in which a fraction of the cells undergoes plasticity; for ε > ε max all cells deform plastically. In the transition region, the fraction of plastic cells P p can be computed as where c and d denote the length of the two segments shown in Fig. 1b. The homogenized stress σ is the average of the stresses in all cells The quantity σ yi (r i ) in the latter equation is a uniformly distributed variable in the interval σ min y , E 0 ε . It follows that the expected value of the expression in Eq. (13) can be calculated as For ε < ε min the response is governed by σ = E 0 ε , while for ε > ε max the material collapses at the constant stress σ = σ max Expected value of lower bound. Now assume that all cells are connected in series and subject to a uniform stress σ . The initial response is linear elastic, σ = E 0 ε , however, as the applied stress attains the value σ = σ min y , the material collapses at constant stress and the macroscopic strain is undefined.
Heterogeneity in hardening modulus. Finally, we consider spatial variations of the hardening modulus H across the material domain according to Eq. (3), while Young's modulus and yield stress are taken as uniform, i.e. �σ y = 0 and E = 0 . In this case, all cells will yield at the same strain ε y = σ y0 /E 0 , as sketched in Fig. 1c.
Expected value of upper bound. All cells are imposed a uniform applied strain ε . In the elastic region ( ε ≤ ε y ) the macroscopic response is dictated by σ = E 0 ε , while for ε > ε y , the homogenized stress σ is the average of the stresses in all cells and given by 2ε . Expected value of lower bound. Now assume that all cells are connected in series and experience the same stress σ . For σ < σ y0 the response is elastic; for σ ≥ σ y0 , the global strain is given by the average of the strains in the cells according to Similar to the analysis presented in "Heterogeneity in Young's modulus" section, the expected value of Eq. (17) can be calculated and the homogenised stress-strain response may be written as which defines the effective strain hardening modulus as �H/ln(H max /H min ).

Numerical calculations
The Monte Carlo Simulation (MCS) method was used to examine the response of elastic-plastic materials with random spatial variations of material properties under different loading conditions. The simulations were carried out in ABAQUS/Standard, conducting repeated FE simulations on 10 realisations of SVEs and evaluating the mean and spread of the populations of outputs. The details of the numerical implementation are described below. FE scheme. Cuboidal or prismatic SVEs of volume V were generated in ABAQUS/CAE via Python scripts and meshed by fully-integrated 8-noded brick elements of cuboidal shape (C3D8) and volume V e . Note that geometric nonlinearity arising from large deformation were accounted for in the analyses. We define as N FE the total number of finite elements in the model N FE = V /V e . A regular tessellation was superimposed to the FE mesh to subdivide the domain into a number N CELL of cuboidal SVEs of volume V CELL , such that N CELL = V /V CELL ; each cell was assigned different material properties to introduce spatial inhomogeneity. The statistical tessellation was such that each cell included an integer number of equally-sized finite elements (hence N CELL ≤ N FE ).
The elastic-plastic microscopic response of the material in each cell was modelled using linear isotropic elasticity (Young's modulus E, Poisson's ratio ν) and incompressible J 2 plasticity with isotropic linear strain hardening (yield strength σ y , hardening modulus H). A pseudo-random number generator in Python was used to generate the values of the microscopic material properties, in accordance with Eqs. (1)-(3).

Boundary conditions and load cases. Uniaxial loading in tension and compression.
We performed simulations of the response in uniaxial tension and compression of prismatic materials domains, which were regularly tessellated and modelled as a set of SVEs. These simulations were performed for comparison to the analytical bounds derived in "Analytical predictions of the uniaxial response" section and to perform numerical convergence analyses. The normal displacements on a face perpendicular to the longitudinal axis of the specimen were constrained, while the opposite face was subject to normal displacements linearly increasing with simulation time, such to induce macroscopic tensile or compressive axial straining. The nodes on the lateral faces were free. Assuming uniform deformation and volume conservation during plastic deformation, the true macroscopic stresses and strains were obtained as and respectively, where the indices 't' and 'n' denote true and nominal stresses or strains, respectively.
Multi-axial loading. We also simulated multi-axial loading conditions to study the effect of strain triaxiality on the response of heterogeneous cubic samples. We imposed uniform displacement boundary conditions on the faces of the cubic domain, forcing these to remain planar as deformation proceeded; as in 17 , the normal displacements of opposite lateral faces were forced to be equal and opposite; these constraints were imposed via the constraint equation tool of Abaqus, making use of appropriate auxiliary nodes. The choice of uniform over periodic boundary conditions was driven by the ease of implementation; this is expected to lead to larger size of RVEs (compared to the case of periodic boundary conditions), however this was not a problem for the simulations in this study, which involved a relatively small number of finite elements. www.nature.com/scientificreports/ In this study we analysed volume elements that were initially statistically isotropic and the applied strains were kept small, such that the anisotropy induced by straining was negligible, as confirmed in preliminary checks. We therefore assumed that the material response could be evaluated in principal strain space, which coincided with principal stress space for our isotropic material. In the FE simulations we therefore prescribed time histories of only the normal macroscopic strains ε xx , ε yy , ε zz , while the macroscopic shear strains were set to zero, ε xy , ε xz , ε yz = 0 (xyz was a reference system aligned with the edges of the cubic domain). In other words, the applied strains ε xx , ε yy , ε zz , were interpreted as principal strains ε I , ε II , ε III . We checked that the macroscopic shear tractions were negligibly small in the simulations, τ xy , τ xz , τ yz ≈ 0 , confirming that the global reference system xyz was a principal system also for the macroscopic stress tensor.
We imposed different values of strain triaxiality by subjecting the cubic domains to normal macroscopic principal strains given by where t represents the simulation time and ε 0 denotes a reference strain rate, which was set to ε 0 = ±0.035 s −1 in all simulations. Nine pairs of parameters (α, β) were chosen, summarised in Table 1, to give 18 different values of strain triaxiality (defined as the ratio of the dilation to the equivalent von Mises strain) in the elastic regime, of which 9 positive ( ε 0 = 0.035 s −1 , loading cases 1-9) and 9 negative ( ε 0 = −0.035 s −1 , loading cases 10-18). Histories of macroscopic true principal stresses and strains were extracted from each simulation from the degrees of freedom and the reaction forces of the auxiliary nodes, and used to determine the corresponding histories of hydrostatic stress σ H and equivalent von Mises stress σ VM , as well as the volumetric strain ε V and equivalent von Mises strain ε VM , as Sensitivity of the uniaxial response to N CELL and N FE . The predictions are affected by the density of the statistical and FE meshes used, representing, respectively, a material length scale and the accuracy of the numerical discretisation. In this section we determine the sensitivity of the predictions to these parameters, with the secondary objective of determining the minimum values of N FE and N CELL which make the predictions approximately insensitive to these same parameters. Note that the results presented in this section were obtained using the loading and boundary conditions described above. Figure  First, we examine the sensitivity of the predictions to N CELL , for the case N FE = N CELL (i.e. each statistical cell is meshed by a single finite element) and a sample with dimensions of 5 × 5 × 5 mm. Figure 2a shows examples of domains with N CELL varying in the range 125-8000. For each choice of N CELL , MCS were performed on two sets of realisations of the microstructure. Set I was generated by implementing spatial variations of the Young's modulus according to Eq. (1), while keeping yield strength σ y0 and hardening modulus H 0 = 0 uniform. For Set II, random variations of the yield stress were considered, according to Eq. (2), keeping elastic modulus and hardening modulus (H 0 = 0) uniform. We define, for Set I, a non-dimensional measure of heterogeneity r E , as the relative variance of the modulus where ΔE and E 0 denote range and average of the random spatial variation of Young's modulus. Such parameter was set to r E = 2 (which represents an extreme case, i.e. the theoretical maximum for r E ) in this set of simulations. Using a similar notation, we define, for Set II, a non-dimensional measure of heterogeneity as A high value r σ = 1.33 was chosen in the simulations. The domains were subject to loading in uniaxial tension, which was interrupted at a total tensile strain of 0.05; the effective macroscopic modulus E eff and the 0.2% proof stress σ eff were extracted from the macroscopic true stress-strain histories, for both Set I and Set II, respectively.  Fig. 2b we summarize the results of these calculations for both Set I and Set II, showing predictions of the normalised effective modulus E eff /E 0 for Set I, and the normalized effective strength σ eff /σ y0 for Set II, for each choice of N CELL . The error bars included in Fig. 2b represent the ranges of each output population (the distributions of outputs were found to be approximately Gaussian, as expected).
The analysis of Fig. 2b shows that σ eff /σ y0 is scarcely sensitive to N CELL , but lower than the value of 1 expected in a "deterministic" simulation (in which material properties are uniformly set at their average values). The spread of such quantity decreases rapidly with increasing N CELL . For Set I, the quantity E eff /E 0 initially increases with N CELL but this dependence becomes small for N CELL > 1000 . Again, the scatter in this normalised measure of strength decreases rapidly with increasing N CELL . Increasing the value of N CELL , while the size of the domain simulated is kept constant, is physically equivalent to increasing the size of the component analysed. Therefore our simulations suggest a size dependence of the tensile response for Set I, such that larger samples are stiffer than smaller ones, as reported by several authors (e.g. 18 ).
Next, we examine the sensitivity of the FE predictions to N FE , with N CELL held constant. FE simulations were performed on a single realisation of the heterogeneous solid, with N CELL = 125 , and the FE mesh was progressively refined, as illustrated in Fig. 3a. Again, two sets of simulations were performed, one with r E = 2 (Set I) and the other with r σ = 1.33 (Set II). Each domain was subjected to uniaxial tension and the simulations were interrupted when a total tensile strain of 0.05 was reached. Figure 3b presents, on the same graph, plots of the normalised effective modulus E eff /E 0 for Set I ( r E = 2 ) and the normalized effective strength σ eff /σ y0 for Set II ( r σ = 1.33 ), for each choice of N FE . It can be seen that both E eff /E 0 and σ eff /σ y0 decrease monotonically with increasing N FE ; the decrease in stiffness is expected, as the number of degrees of freedom of the calculation increases; the decrease in strength is justified by the fact that a progressively refined FE mesh better captures the stress concentrations induced by the microscopic variations of mechanical properties.  www.nature.com/scientificreports/ The latter mesh sensitivity study was repeated for higher values of the statistical mesh density, namely N CELL = 1000 and N CELL = 8000 , for both Set I and Set II. The predictions are presented in Fig. 4, where E eff /E 0 is plotted as a function of the ratio N FE /N CELL in Fig. 4a for Set I ( r E = 2 ), while for Set II ( r σ = 1.33 ), predictions of σ eff /σ y0 are shown Fig. 4b

Results and discussion
Uniaxial response. In this section we analyse the numerical predictions of the uniaxial response in detail and compare them to the bounds derived analytically above. We chose prismatic volume elements (VEs) of square cross-section, consisting of 20,736 cubic cells (N CELL = 20,736) and each cell was meshed with a single finite element ( N FE = N CELL ). The prismatic domains analysed had dimensions of 6 × 4 × 4 mm , to represent the gauge portion of test specimens typically used for laboratory-scale mechanical tests. As shown in Fig. 2b, this choice of N FE , N CELL makes the predictions of average modulus and strength insensitive to further refinements of the discretisation. All the numerical predictions presented in this section are average responses obtained from 10 simulations of different realisations of the heterogeneous specimens.
Effect of variation of Young's modulus. We now present the results of simulations in which a random spatial variation of Young's modulus was imposed, with uniform yield strength and uniform hardening modulus, H = 0 (as sketched in Fig. 1a). The average macroscopic uniaxial true stress versus true strain histories are presented in Fig. 5a for the case of tension and in Fig. 5b for the case of compression. Simulations were conducted at three different values of r E as indicated; the insets show contours of the absolute value of the maximum principal strain at a 5% total macroscopic axial strain (for the maximum value of r E considered). A "deterministic" simulation, in which all material properties were taken as uniform and equal to the average values (E 0 = 70 GPa; σ y0 = 300 MPa; H = 0) is included in the figures for comparison.
In both tension and compression, a non-linear elastic-plastic response begins at macroscopic stress well below σ y0 = 300 MPa , as a consequence of the stress concentrations and multiaxial stress states induced by the variations of Young's modulus; the stress at the onset of plasticity decreases with increasing r E . At higher stresses, the material exhibits a non-linear macroscopic response reminiscent of strain-hardening, but non-linearity is rather a consequence of the progressive yield of more and more cells, similar to what described by Eq. (8) (we recall that the microscopic material response is perfectly plastic). The peak true stress is scarcely sensitive to At relatively large strains, the mechanisms of deformation are different in tension and compression. In tension the deformation, initially uniform, localises and a neck forms in the sample at a random location, as evident from the inset of Fig. 5a. Due to the local reduction of cross-sectional area in the neck, the macroscopic response exhibits a stress peak followed by geometric softening. Furthermore, the macroscopic strain at which a peak in tensile stress is attained is sensitive to r E : the higher r E , the higher the strain at peak stress; this suggests that the ductility of a heterogeneous material increases with the degree of elastic heterogeneity, as quantified in this case by the parameter r E . In compression, on the other hand, the necking mechanism is absent and the macroscopic response progressively evolves towards a plateau response at stress σ y0 = 300 MPa , as shown in Fig. 5b. Some inhomogeneity in the strain field is visible, due to the tendency of the specimen to localise deformation along inclined planes. In the apparent hardening phase of both the tensile and compressive responses, the flow stress at any given macroscopic strain decreases with increasing r E . Figure 5c presents a comparison between the numerically predicted normalised effective tensile modulus and the expected value of the analytical bounds developed above. It can be seen that the average values of our numerical predictions fall between the analytically predicted expected values of the upper and lower bounds.
Effect of variation of yield stress. In Fig. 6a,b we present average macroscopic true stress versus true strain histories obtained from simulations in which the yield stress was varying in space, while the Young's modulus and hardening modulus were kept uniform E 0 = 70 GPa; σ y0 = 300 MPa; H = 0 . Predictions are presented for three different levels of heterogeneity, quantified by the values of r σ indicated. The insets show contours of the absolute value of the maximum principal strain at 5% total macroscopic axial strain.
The effects of a microscopic variation of yield stress are somewhat different from those observed for variations in Young's modulus. In both tension and compression, compared to the deterministic simulation, the macroscopic response of the heterogeneous material features an onset of plasticity at relatively low values of stress, followed by a non-linear (apparent) strain hardening phase, consequence of the progressive yielding of the cells in the sample. Again, a higher degree of heterogeneity r σ results in lower flow stress at any applied macroscopic strain. In tension and compression, deformation localises along shear bands at approximately ±45° on the loading axis; in tension these bands evolve to trigger a macroscopic necking mechanism associated with strain softening, while in compression, the predicted stress quickly reaches a plateau after the initial elastic phase, as shear bands spread across the whole sample. In tension, the macroscopic strains at necking increase with r σ , reinforcing the notion that material heterogeneity increases macroscopic ductility.
In Fig. 6c numerical predictions of the 0.2% proof stress in tension are compared to the analytical bounds. As expected, these lie within the analytical bounds but appear much closer to the upper bound than to the lower bound.
Effect of probability distribution. We performed additional MCSs using a symmetric, triangular (T) probability density function of yield stress, to compare with our previous results obtained for the case of a uniform, rectangular distribution (R) with r σ = 1.33 (Fig. 6). The range of variation for the triangular probability density, r ′ σ = �σ ′ y /σ ′ y0 , was chosen such to obtain the same average and variance as for the uniform distribution with r σ = 1.33 . This was achieved by setting r ′ σ = √ 2r σ = 1.88 and σ ′ y0 = σ y0 = 300 MPa . In Fig. 7 we compare the average macroscopic stress-strain responses obtained with the uniform distribution ( r σ = 1.33 ) to those obtained with a triangular distribution of equivalent variance ( r ′ σ = 1.88 ). Clearly the two types of predictions www.nature.com/scientificreports/ are nearly coincident in both tension (Fig. 7a) and compression (Fig. 7b), suggesting a scarce sensitivity of our predictions to the assumed shape of the probability density functions which govern the spatial variation of material properties. It would be interesting to examine the case of a non-symmetric probability density function of the same average and variance, which is omitted here for brevity.

Effect of variation of hardening modulus.
A separate set of simulations was conducted by imposing a spatial variation of the material hardening modulus, with elastic modulus and yield stress kept constant E 0 = 70 GPa; σ y0 = 300 MPa . The mean value of the hardening modulus was H 0 = 35 Gpa and the degree of variation was chosen as r H = 2.
The ensuing macroscopic response, not shown for the sake of brevity, comprised a linear elastic part followed by a linear strain-hardening phase. The effective hardening modulus, for the (relatively high) degree of variation r H considered, was much closer to that predicted by the upper bound developed above (Eq. (16)) than to the corresponding lower bound (Eq. (18)), indicating cooperative straining. A spatial variation of the hardening modulus alone does not have substantial effects on the macroscopic response, with the exception of resulting in a macroscopic hardening modulus only slightly smaller than the average microscopic value (as predicted by the expected value of upper bound, Eq. (16)).
Size (and shape) dependence of the response. Introducing a spatial variation of material properties in the numerical models introduces a length scale in the problem, namely the dimension of the SVE used. Dimensional analysis dictates that the problem must depend on a non-dimensional parameter involving this cell dimension, for example the ratio between size of the material sample and size of the statistical cells, which increases mono-  www.nature.com/scientificreports/ tonically with decreasing N CELL . The number of SVEs (or cells) used, N CELL , is representative of physical size of the material specimen. The dependence of the macroscopic response upon N CELL has been already highlighted in Fig. 2. We now examine in detail how the macroscopic stress-strain response changes as a function of N CELL . First, consider the case of uniform elastic modulus, uniform hardening modulus and non-uniform yield stress (r σ = 1.33 , r E = r H = 0) . Figure 8a,b present, respectively, the tensile and compressive stress-strain curves for a heterogeneous material, compared to the deterministic (homogeneous) case. Each figure includes three curves corresponding to different values of N CELL . We include the corresponding value of the ratio ϕ , defined as the number of material cells lying on the specimen's lateral surface divided by those lying within the specimens. The latter ratio depends on N CELL through a simple geometric relation and it also depends on the shape of the specimen analysed. It is seen from Fig. 8 that the material displays an inverse size effect in both tension and compression, with larger specimens (or equivalently, specimens with larger values of N CELL ) displaying a stronger material response. This has been reported by other studies on size effects in the response of ductile elastic-plastic solids (e.g. 19 ) and is justified in terms of the different stress states of material domains lying on the surface of the solids and those contained within this surface. In fact, material cells lying with at least one side on the external surface are less constrained, having at least one vanishing principal stress, and tend to yield earlier than material cells within the bulk of the SVE. Our results support this notion, reporting a weaker response in both tension and compression for higher values of ϕ . We note that assuming a uniform material response would not give rise to such size effect, but modelling the heterogeneity in material properties results in a dependence of the macroscopic response upon both size and shape of the specimen analysed.
In the data presented in Fig. 8 the aspect ratios of the specimens analysed were kept constant. In a separate set of simulations (r σ = 1.33 , r E = r H = 0, N FE = N CELL ) , we varied the geometry of the SVE substantially, ranging from a chain of cells in series to an array of cells all arranged in parallel; the corresponding uniaxial tensile strength varied from the mean yield stress (for the case of a parallel arrangement, in line with the upper bound derived in "Analytical predictions of the uniaxial response" section) to approximately one third of this value, for  www.nature.com/scientificreports/ the case of a chain arrangement. Such variations are driven by the dependence of ϕ upon the geometry (aspect ratios) of the SVE, but also by different mechanisms of deformation and different stress triaxiality induced by the specimen geometries. Next we present, in Fig. 9, analogous information as that in Fig. 8, for the case of a spatial variation of the hardening modulus alone (r H = 2 , r E = r σ = 0) ; the Young's modulus and yield strength were chosen as E 0 = 70 GPa and σ y0 = 300 MPa , respectively, and the mean value of hardening modulus was H 0 = 35 GPa . The observed size effect is in this case very mild in the wide range of sizes ( 324 ≤ N CELL ≤ 20,736 ) considered here.

Multiaxial response.
To illustrate the effects of material heterogeneity on the multiaxial response, we used the loading and boundary conditions described in the "Multi-axial loading" section to construct pseudo-yield surfaces of heterogeneous solids and compare with the homogeneous case. Cubic SVEs were generated with N FE = N CELL = 21,952 , using a 5 × 5 × 5 mm specimen. Macroscopic strain histories were imposed on the SVEs according to Eq. (19) and Table 1. Two sets of simulations were conducted, varying either the elastic modulus or the yield stress at microscopic scale. All pseudo-yield surfaces were constructed using the average responses obtained from 10 simulations of different realisations of the heterogeneous cubic domain. Figure 10 shows the yield surface of a heterogeneous material featuring microscopic variations of the elastic modulus. In Fig. 10a we present, for r E = 1.43 , the evolution of the yield surface. Trajectories of the hydrostatic and deviatoric stresses at different (initial) stress triaxiality are shown, and on each trajectory two points are marked, corresponding to macroscopic von Mises equivalent strains of 0.01 (squares) and 0.04 (circles). Clearly the macroscopic yield surface is pressure-dependent, although the microscopic material response is pressure-insensitive (obeying the von Mises criterion); this is due to the fact that an imposed macroscopic hydrostatic strain results in deviatoric stress components at microscopic level, due to the elastic heterogeneity. The macroscopic material hardening can be deduced from the figure and is neither isotropic nor kinematic (note that isotropic hardening is imposed at microscopic level). Figure 10b shows the effect of the degree of heterogeneity r E on the multiaxial response: clearly, a higher heterogeneity gives a weaker response and enhances the pressure-sensitivity of the macroscopic response (note that the yield surfaces in Fig. 10b are obtained at macroscopic von Mises equivalent strains of 0.01). Figure 11 presents similar information for the case of microscopic variations of the yield stress, with uniform elastic modulus; Fig. 11a shows the evolution of the yield surface for r σ = 1.33 , while Fig. 11b presents the effect of varying r σ on the macroscopic yield surface. In this case the macroscopic response is pressure-insensitive and hardening is isotropic, as for the local constitutive response; this is due to the fact that the material is elastically homogeneous and macroscopic hydrostatic stress components do not induce microscopic deviatoric components in this case. Again, a higher degree of heterogeneity results in a weaker macroscopic response (Fig. 11b).
Can we predict the fatigue endurance limit from the quasi-static response? The discussion of the results presented above has highlighted, among other things, the fact that in real materials the onset of yield occurs at relatively low stresses (see Figs. 5 and 6), much lower than what is commonly referred to as the material 'yield stress' (i.e. proof stress, or macroscopic flow stress at 0.2% plastic strain). In addition, we found that the macroscopic response of heterogeneous materials featuring variations of the yield stress is close (Fig. 5) to the predictions of the simple upper bound model developed in "Analytical predictions of the uniaxial response" (Eq. (14)). In principle, fitting Eq. (14) to the stress strain curve of a material with small macroscopic strain hardening may be used to determine the macroscopic stress at the onset of local plasticity, σ min y . Such fitting exercise was conducted using uniaxial tension data by Bettaieb et al. 20 on two types of Titanium alloys, considering only the stress-strain histories for strains less than 0.01; an initial linear response was assumed, followed by a non-linear response (given by Eq. (14)) for σ > σ min y . The value of σ max y was taken as the flow stress at total strain of 0.01 (this value approximately corresponds to the point at which the measured stress-strain curves in 20 attain a horizontal tangent). This fitting exercise yielded, for the two alloys considered (denoted as Ti-5553-1 and Ti-5553-3 in 20 ) σ min y = 334 MPa and σ min y = 312 MPa, respectively. The authors also measured S-N curves for these two alloys (at relatively low stress ratio, R = 0.1) and found that both alloys www.nature.com/scientificreports/ displayed a fatigue endurance limit, measured as 330 MPa and 310 MPa for Ti-5553-1 and Ti-5553-3, respectively. These are remarkably close to the values of σ min y extracted from the uniaxial stress-strain curves. It is well known that fatigue is a macroscopic phenomenon initiated by microscopic plasticity at relatively low levels of stress, and due to local deformation and fracture mechanisms supported and enhanced by cyclic loading. Although the mechanisms are quite complex, and surely they involve complex loss of cohesion, it is logical to think that if the applied stress levels are low enough to not initiate noticeable microscopic plasticity, then the material will not fail by fatigue at such applied stresses. The data, in this case, support this notion; however more complex modelling techniques (accounting for material strain hardening, loss of cohesion, more sophisticated bounds) and a more extensive experimental campaign on different materials would be needed to substantiate this hypothesis.
However, to further support the possibility of a link between the fatigue endurance limit of metals and their quasi-static monotonic response, we analyse data from a material database 21 , as shown in Fig. 12. The database stores mechanical properties for a large number of engineering materials. In particular, for unreinforced metallic alloys, it reports (i) the yield stress σ 0.2 , intended as the stress at a plastic strain of 0.2%, (ii) the tensile strength σ T , intended as the peak stress in uniaxial tension, (iii) the Young's modulus E 0 , and (iv) the fatigue strength at 10 7 cycles σ 10 7 F , which is in practice considered as the fatigue endurance limit. Figure 12a-c show the correlation between such fatigue strength and elastic modulus, yield strength and tensile strength, respectively, for the 186 materials for which the four properties above are measured via mechanical test (rather than estimated, by approaches similar to that in 22 ). The fatigue strength in general increases when these three parameters increase, and in first approximation it does so according to power-laws of exponents close to 1.
We now assume that the stress versus strain response of metallic alloys can be described by Eq. (14), implying a negligible heterogeneity in elastic properties, moderate variability of the yield stress and negligible strain hardening modulus; the proof stress σ ε at a certain plastic strain ε (in this case σ 0.2 , evaluated at plastic strain ε = 0.002 ) can be determined by finding the intersection of Eq. (14) with the line of equation σ = E 0 (ε − ε) . Under our hypothesis, we assume that σ 10 7 F = σ min y . With regards to σ max y , we assume that for materials with small strain hardening modulus it is σ T ≈ σ max y ; the logic of this is that for a material with a locally perfectly-plastic www.nature.com/scientificreports/ response, global plastic collapse (and consequent ultimate failure) will have occurred when the applied stress attains the yield stress of the strongest SVE. This leads to the equation Equation (28) makes mathematical sense for σ T ≥ σ 10 7 F , which is satisfied for all materials in the dataset investigated. To further ensure that it can be solved for σ 10 7 F , it must also satisfy (σ T − σ 0.2 )/E 0 ≥ ε/2 . In Fig. 12d we plot measured data for a subset of 149 metallic alloys satisfying the condition above. It is clear that the quantities at the two sides of Eq. (28) are correlated. In Fig. 12e we present the same information as in Fig. 12d, but imposing the further condition σ T < 1.2σ 0.2 , to exclude materials with very pronounced strain hardening. The correlation between the two sides of Eq. (28), including this time only 60 materials, is plotted in linear scale (rather than logarithmic), to highlight the discrepancies between data and prediction (included in both Fig. 12d,e). Again, Eq. (28) provides a good description of the measured data.
In consideration of the enormous impact that correlations similar to Eq. (28) would have in engineering practice, we think that the preliminary findings highlighted in this section should be reinforced by more extensive investigations, and we hope that some of the readers will join us in this.

Concluding remarks
Monte Carlo analyses and FE simulations were used to study the effects of microscopic variations of elastic modulus, yield strength and hardening modulus upon the macroscopic elastic-plastic response of model heterogeneous materials. Volume elements consisting of regular arrays of SVEs were generated and their uniaxial and multiaxial stress-strain responses were predicted. The mechanical properties of the SVEs were taken as linear elastic, followed by incompressible J2 plasticity with constant strain hardening modulus and isotropic hardening. Mechanical properties were different in each SVE, according to uncorrelated random fields of uniform probability density. Microscopic variations of either the elastic modulus, the yield stress or the hardening modulus were studied, and their macroscopic effects were investigated. The main conclusions and observations of the study are as follows: • The macroscopic elastic modulus and the yield stress extracted from experiments are lower than the spatial averages of their microscopic counterparts. • Heterogeneous materials with a local perfectly-plastic response may display an apparent macroscopic strainhardening. • The macroscopic response of a heterogeneous solid is scarcely sensitive to the shape of the probability density function governing the microscopic spatial variability of the mechanical properties, but depends strongly on the relative variance of such probability density. • The onset of yield in a heterogeneous material occurs at stresses much lower than the proof stress; the higher the degree of heterogeneity, the lower is the macroscopic stress to cause the first local yield. • Microscopic variation of stiffness or yield stress induce tension/compression asymmetry, strain localisation by shear banding and necking instability in uniaxial tension. • The strain at which tensile necking occurs in ductile materials increases with the degree of heterogeneity.
• Heterogeneity in plastic properties results in an inverse size effect as well as in a pronounced dependence of the response on the shape of the specimen. • An elastically heterogeneous material, comprising an array of plastically incompressible domains with isotropic hardening, displays a pressure-sensitive macroscopic response with non-isotropic hardening. • A preliminary investigation on the application of the findings of this paper suggests that the fatigue endurance of solids might be estimated from a detailed analysis of the stress versus strain curves measured in quasi-static monotonic experiments.