Fracture initiation in multiphase materials: a systematic threedimensional approach using a FFTbased solver
Abstract
This paper studies a twophase material with a microstructure composed of a hard brittle reinforcement phase embedded in a soft ductile matrix. It addresses the full threedimensional nature of the microstructure and macroscopic deformation. A large ensemble of periodic microstructures is used, whereby the individual grains of the two phases are modeled using equisized cubes. A particular solution strategy relying on the Fast Fourier Transform is adopted, which has a high computational efficiency both in terms of speed and memory footprint, thus enabling a statistically meaningful analysis. This solution method naturally accompanies the regular microstructural model, as the Fast Fourier Transform relies on a regular grid.
Using the many considered microstructures as an ensemble, the average arrangement of phases around fracture initiation sites is objectively identified by the correlation between microstructure and fracture initiation – in three dimensions. The results show that fracture initiates where regions of the hard phase are interrupted by bands of the soft phase that are aligned with the direction of maximum shear. In such regions, the hard phase is arranged such that the area of the phase boundary perpendicular to the principal strain direction is maximum, leading to high hydrostatic tensile stresses, while not interrupting the shear bands that form in the soft phase. The local incompatibility that is present around the shear bands is responsible for a high plastic strain. By comparing the response to a twodimensional microstructure it is observed that the response is qualitatively similar (both macroscopically and microscopically). One important difference is that the local strain partitioning between the two phases is overpredicted by the twodimensional microstructure, leading to an overestimation of damage.
Fracture initiation in multiphase materials: a systematic threedimensional approach using a FFTbased solver
Keywords: micromechanics; ductile failure; damage; multiphase materials; FFT solver
1 Introduction
Multiphase materials are used in a wide variety of applications because they often encompass good engineering properties: electrical, thermal, mechanical, etc. By combining two or more distinct phases at the level of the microstructure, a combination of properties is obtained that cannot always be achieved otherwise. This paper focuses on the class of twophase materials that display both strength and ductility using a microstructure in which a soft and ductile matrix is reinforced using a hard and brittle phase.
The elastoplastic response of such materials is reasonably well understood. Simple rules of mixture yield satisfactory predictions for many purposes [1, 2], but also meanfield [3, 4] and unit cell models have been applied [5, 6, 7, 8]. Many of them have been extended to incorporate microstructural damage to predict the loss of material integrity observed at the macroscale. However, in most cases, the role of the microstructure, and in particular of the phase distribution, is only partially considered and not fully understood.
Different researchers have performed studies that aimed at characterizing the role of the spatial distribution of the phases on the initiation and growth of fracture at the level of individual grains [9, 10, 11, 12]. In these papers it was found that damage preferentially nucleated in or inbetween reinforcement particles that are neighboring along the tensile axis – see also the experimental work by Lewandowski et al. [13]. A positive hydrostatic stress was one of the configurational characteristics driving damage nucleation and growth.
Recently, De Geus et al. [12] have demonstrated that the role of the microstructure can be rationalized by considering the average phase distribution around damage sites through a conditional spatial average. The strength of this approach is that the results need little interpretation as the correlation between morphology and fracture initiation is quantified. Besides confirming the role of neighboring reinforcement particles in the tensile direction, the critical microstructure features a pattern composed of matrix phase directly adjacent to the damage site in the direction perpendicular to the tensile direction, and in bands in the directions of maximum shear. Later, it was shown that this phase distribution is also critical for fracture initiation by cleavage in the reinforcement phase [14].
Most published models are restricted to two spatial dimensions. Likewise experimental observations generally are restricted to crosssections. In spite of the valuable insights gained, the impact of the geometrical approximation remains unclear. Some studies were performed in three dimensions [10, 11, 15], but a systematic analysis including many statistical variations is often unfeasible. Ramazani et al. [15] observed a clear difference in the macroscopic mechanical properties predicted by two and threedimensional micromechanical models. For a different physical process, Xu et al. [16] compared the accuracy of twodimensional representative volume elements (RVEs) to full threedimensional RVEs in the context of humidity induced electrochemical degradation of fuel cells with a Nickel – Yttriastabilized Zirconia microstructure. They concluded that the yield stresses are of the same order of magnitude in the twodimensional model as in the threedimensional model, but the former cannot capture the thermoelectrical properties accurately as its percolation is inherently limited. These distinct mechanisms may also have a strong influence on fracture initiation for the class of materials considered here, as the plastic strain localizes in a path of soft phase that is connected in the direction of shear, prior to fracture initiation.
A systematic analysis of the local mechanical response of a threedimensional microstructure, with statistical relevance, is still beyond the capabilities of conventional computational solid mechanics, based on the finite elements method. Different methods have been developed that enable such an analysis. Among them, the Voronoi tesselation finite elements method of Moorthy and Ghosh [17] is worth mentioning. This method employs a microstructure that is generated using a Voronoi tessellation. Each Voronoi cell forms an element for which the response is calculated using an assumed stress distribution. Although the resulting number of degrees of freedom is limited, the extension to three dimensions is still nontrivial [18].
Besides the different strategies for improving the efficiency of the calculations based on the finite elements method, fullfield calculations based on the Fast Fourier Transform (FFT), as first proposed by Moulinec and Suquet [19], have proved to be a serious alternative for periodic inhomogeneous media when the number of degrees of freedom becomes significant, as it is the case in 3D. This method takes advantage of the fact that equivalent linear elastic problems can be formulated and recast into integral forms, where the relevant Green functions are convolved with some appropriate source fields. The convolution product is expressed as a simple product in Fourier space and is computationally evaluated by means of an efficient FFTlibrary. Even though the application of the Fast Fourier Transform restricts the discretization to a regular grid, this method has the great advantage that the computational costs, both in terms of CPU time and memory footprint, scale approximately linearly with the number of grid points. The method can thus handle large systems and appears to be the best suited for 3D microstructures.
The current paper aims to extend the analysis by De Geus et al. [12] to a threedimensional setting. By comparing different threedimensional deformation paths for a large number of threedimensional microstructures, the role of the microstructural morphology on fracture initiation is systematically identified. The results are compared to those obtained from twodimensional microstructures, complementing the analyses of fracture initiation mechanisms but also allowing to identify the limitations of such twodimensional models. This is achieved by means of the algorithm proposed by Moulinec and Suquet [19] to solve mechanical equilibrium in heterogeneous materials. Besides its low memory footprint, the simplifying but reasonable assumption that the elastic moduli are homogeneous renders the method incomparably fast and the most suitable for our statistical analysis.
This paper is structured as follows. The microstructural model and the numerical implementation are discussed in Sections 2 and 3 respectively. In Section 4, the two and threedimensional microstructures are compared in terms of macroscopic and microscopic elastoplastic response and fracture initiation. For the threedimensional microstructures, the effects of different deformation paths are investigated in Section 5. The paper ends with several concluding remarks in Section 6.
Nomenclature
second order tensor  
fourth order tensor  
dyadic tensor product ()  
single tensor contraction ()  
double tensor contraction ()  
ensemble average  
volume (unit cell) average 
2 Model
2.1 Microstructure
The microstructure is represented by an ensemble of threedimensional periodic unit cells. Each unit cell comprises cubes arranged in a regular grid, see for a typical example Figure 1(b). These cubes represent the individual grains or particles of the two phases considered – a soft and a hard phase with respectively low and high yield stresses. They are referred to as grains below, to avoid confusion. It must be emphasized that the model does not account for any crystalline orientation such that the boundaries between grains of the same phase are mechanically meaningless. The introduced idealization enables a clearly defined statistical analysis in which the average spatial distribution of phases around a fracture initiation site can be transparently quantified. It must be stressed that to provide a good estimate of the strain and stress fields, the numerical discretization is finer than the grid of grains, as discussed in Section 3. However, consistent with the adopted idealization, only grain averaged quantities are considered.
The phases are randomly distributed in each unit cell according to a predetermined volume fraction of the hard phase. To achieve this, each individual grain in each of the unit cells is assigned the properties of the hard phase if a random number in the interval is smaller than the chosen . All other grains are assigned the properties of the soft phase. Throughout this paper, is used. The obtained ensemble average volume fraction of the hard phase, , is within of the target value . The hard phase volume fraction in the individual unit cells varies between around this value.
2.2 Constitutive model
Different constitutive models provide a satisfactory description of the micromechanical response at the considered scale and level of approximation, that of an aggregate of grains. Here, both phases are modeled using small strain elastoviscoplasticity, whereby the parameters for the viscous part are taken to effectively approximate a rateindependent model. It has been verified that this model is sufficiently accurate for the purpose of this study by comparing the response of the selected model to rateindependent models in small and finite strain; the results of this comparison are not included here for conciseness.
In the constitutive model, the linear strain tensor is additively split in an elastic and a plastic part:
(1) 
with the elastic strain tensor and the plastic strain tensor.
Assuming linear elasticity, the stress tensor is determined from the elastic strain tensor by Hooke’s law:
(2) 
with the Cauchy stress tensor. is the elastic stiffness, which is assumed to be homogeneous (independent of space) and isotropic:
(3) 
where the bulk modulus and the shear modulus depend on the Young’s modulus and Poisson’s ratio in the usual way: and . The tensors and respectively are the secondorder identity tensor and the fourthorder symmetric identity tensor.
Plasticity is modeled using standard plasticity, by adopting the following yield function:
(4) 
with Von Mises’s equivalent stress – defined as , with the stress deviator – and the yield stress. The latter depends on the accumulated equivalent plastic strain
(5) 
using a linear hardening relation:
(6) 
with the initial yield stress and the hardening modulus.
Associative flow is assumed to model the evolution of the plastic strain rate:
(7) 
where the plastic multiplier is obtained through Norton’s law, i.e.
(8) 
with the rate dependence exponent and a material constant.
The two phases have the same elastic moduli and differ only through their plastic behavior: the hard phase has a yield stress and hardening modulus which are twice those of the soft phase. The material parameters are summarized in Table 1. Note that the ratedependence exponent of adequately approximates a rateindependent elastoplastic model.
2.3 Applied deformation
As a result of the assumed periodicity of the microstructure in all spatial directions, only the volume averaged deformation needs to be prescribed. In this way, local (microscopic) fluctuations in stress and strain are permitted throughout the cell and along its boundary. In the following, two different average strain paths are considered: planar and axisymmetric shear. Both strain paths are isochoric and result in a vanishing macroscopic stress triaxiality. Note that the triaxiality in the individual grains is determined by equilibrium and compatibility, and generally does not vanish.
The applied macroscopic strain tensor, , for planar shear reads
(9) 
in which is the macroscopic equivalent strain. Note that all directions are treated in the same way, fluctuations are thus allowed in the direction equivalently to the  and direction. Locally the deformation is therefore threedimensional.
For axisymmetric shear, the strain tensor reads
(10) 
where retains its meaning of the macroscopic equivalent strain.
For both cases, the microstructure is deformed up to a value of . This value is large enough to be in a relevant regime for fracture initiation, but small enough to stay within the limits of the small strain approximation.
2.4 Fracture initiation criterion
The initiation of fracture at the local grain level is modeled using a damage descriptor. This descriptor reveals fracture initiation in the grains, but does not cause any material degradation. For the soft phase the JohnsonCook model is used, which is commonly applied to study ductile fracture initiation [20]. In this model the rate of accumulated plastic strain is compared to a critical plastic strain and then integrated over the deformation history, i.e.
(11) 
The accumulated plastic strain rate, , is the time derivative of (5) and is therefore nonnegative; the critical plastic strain is a function of the (local) stress triaxiality as follows:
(12) 
The the stress triaxiality is the ratio of the hydrostatic stress and the equivalent stress ; the material parameters , , and are summarized in Table 1. Fracture initiation occurs in the grain once the descriptor . It is assumed that there is no fracture initiation in the hard phase.
2.5 Fracture initiation hotspot
The relation between the microstructure and initiation of fracture is studied by calculating the characteristic phase distribution around the fracture initiation sites. To obtain a statistically representative result, all fracture initiation sites in the entire ensemble of the random unit cells are taken into account at once. The resulting average pattern indicates what is the probability to find hard and soft phase in the vicinity of any fracture initiation site. This method, as introduced in [12], is discussed in detail based on a single unit cell. The extension to the ensemble average is straightforward and therefore omitted.
A phase indicator, , is introduced to define the phase for each grain as a threedimensional scalar field, i.e.
(13) 
where corresponds to the position in the unit cell. For the cubic microstructures considered in this paper, the index is the voxel index of each grain. Similarly, fracture initiation is indicated by which is unitary if the damage descriptor , and equals zero in all other cases.
The (spatially) average(d) phase distribution around the fracture initiation sites (where ) is given by the following weighted average phase probability function:
(14) 
where is the relative position with respect to the initiation fracture sites, and where in the sums loops over all grains in the unit cell, taking the periodicity into account. Note that (14) effectively quantifies the crosscorrelation between and , and that it may be evaluated as a convolution (e.g. with FFT). It may be interpreted as the probability of finding the hard phase as a function of the position relative to the fracture initiation site. As the microstructure consists of only two phases, the probability of finding the soft phase is also contained in this result: if at a certain relative position (the hard phase volume fraction), then hard phase at this relative position promotes fracture initiation, while if , soft phase at that relative position will promote fracture initiation.
3 Spectral solver for mechanical equilibrium
3.1 Introduction
The numerical scheme introduced by Moulinec and Suquet [19] is used in this paper. In this scheme, mechanical equilibrium is solved iteratively whereby a reference elastic medium is used to project an intermediate solution to a new, compatible, intermediate solution until equilibrium is satisfied. This method does not require the assembly of a tangent matrix or the solution of a linear system, and therefore its memory footprint scales approximately with the number of gridpoints instead of its square. The projection itself is performed using the Discrete Fourier Transform which enhances the computational efficiency. Although it has been shown that the convergence rate of this particular algorithm depends critically on the choice of the reference medium, e.g. [24], this issue is avoided in the present study because the microstructure is considered as elastically homogeneous. For the same reason, oscillations in the solution (e.g. the wellknown Gibbs phenomenon) are also avoided here.
Below, we recall briefly the algorithm of Moulinec and Suquet [19], before explaining how the specific constitutive model enters this algorithm.
3.2 Solution using auxiliary problem
For this method, the local strains are the unknowns. The goal is therefore to find a field of periodic strains, , that satisfies mechanical equilibrium, i.e.
(15) 
An auxiliary problem is introduced for which the stress is the response of a homogeneous linear elastic body, with stiffness , subjected to a polarization stress , i.e.
(16) 
The solution to this problem is given by the following convolution, known as the periodic LippmannSchwinger equation:
(17) 
with the periodic Green operator associated with and the average strain tensor. The convolution can be evaluated in Fourier space using a simple tensor contraction, as follows:
(18) 
where the denotes that quantities have been transformed to Fourier space, and denotes the wavevector dependence in that space.
The polarization stress is generally composed of two contributions: (i) one ensuing from the inhomogeneity of elastic constants and (ii) one involving all eigenstrains such as the plastic strain :
(19) 
Since both and are not known yet, the solution has to be found iteratively. The rate of convergence is thereby determined by the choice of the reference medium (through ).
3.3 Problem specific solution
In the present work, the microstructure is elastically homogeneous, i.e. independent of , and the polarization stress becomes:
(20) 
where denotes the current time increment. Hence, iterations can be entirely avoided by using an explicit time integration scheme for assessing the plastic strain at a new time step:
(21) 
(22)  
(23) 
Finally, knowing the polarization stress, its Fourier transform is inserted into (18) to determine the Fourier transform of the strain field at the new time increment, . The inverse Fourier transform then provides the corresponding strain field in real space, , which were the unknowns.
3.4 Applied discretization
The Fourier transform and its inverse are evaluated using the Discrete Fourier Transform, for which efficient open source implementations are readily available. Each of the grains of the microstructure is discretized using grid points (see Figure 1(a)). It has been verified that this discretization leads to a converged solution in terms of the grain averaged tensors used in this paper. The refinement study confirmed that the approximation is smaller than 1% relative error in stress and (plastic) strain measured in the individual grains.
To minimize the approximation of the explicit time integration, the deformation is applied in small steps. It has been verified that this temporal discretization is small enough to ensure stability and convergence of the scheme. Note that even though many steps are applied, the computations at each time step are fast, even for the large threedimensional unit cells.
4 Twodimensional vs. threedimensional microstructure
The effect of the threedimensional (3D) nature of the considered microstructures is investigated by comparing them to their twodimensional (2D) counterpart for which each crosssection in direction is the same. In this comparison both the 2D and 3D microstructures are subjected to macroscopic planar shear (see Equation 9). Since the number of considered morphologies is smaller for the 2D microstructures, the size of the ensemble is increased to random unit cells.
4.1 Macroscopic response
The computed macroscopic equivalent stress response, , as a function of the applied equivalent strain, , is shown in Figure 2. For the 3D microstructures, in black, and the 2D microstructures, in blue, the upper and lower bounds of the different unit cells are plotted. The constitutive response of the individual phases is also shown using a dark blue curve for the soft matrix and red for the hard reinforcement phase. As observed, all curves coincide in the elastic regime, as the microstructure is elastically homogeneous. In the plastic regime, the response of the microstructure is a nonlinear combination of the response of the individual phases. The spread between the responses of the different unit cells is larger for the 2D model than for the 3D model. The fact that the scatter in volume fraction between the 2D unitcells is somewhat larger (it is in the range ) thereby only partly explains the difference. A clear difference is also observed in the plastic response directly after yielding. The response of the 2D model is well below that of the 3D model. The difference is explained by examining the local plastic response in the next section.
4.2 Microscopic response
The local accumulated plastic strain, , is shown in Figure 3 for the 2D microstructure and in Figure 4 for the 3D microstructure; in both cases for an applied equivalent strain . The figures include a threedimensional view of the unit cell and planar views of the visible (outer) planes (for the 2D microstructure only one view is relevant). Note that because of the periodicity of the microstructure, the outer planes are in fact crosssections of the material.
For both 2D and 3D microstructures, it is observed that the highest values of occur in the soft phase (with an average of and respectively). In most of the hard grains is considerably lower (with an average of and respectively). From the crosssection along the direction (Figures 3(b) and 4(c)), it is observed that the highest values of are localized in bands of connected soft phase at degree angles with respect to the axis, in particular if these bands are surrounded by hard phase. For the 3D microstructure, the crosssections along the plane and the plane (Figure 4(b,d)) reveal that these shear bands are extended along the axis.
Plastic strain thus localizes in shear bands at degree angles with respect to the  and axis. This coincides with the direction of maximum shear of the applied macroscopic strain. It can be concluded that the highest plastic strain occurs where the soft phase is connected in the direction of shear, and surrounded by the hard phase.
These observations, made on a single unit cell, can be extended to the entire ensemble. In Figure 5 the probability density, , of the accumulated plastic strain, , is shown for different strain increments. The figure includes both the 2D (Figure 5(a–b)) and 3D microstructures (Figure 5(c–d)); the hard phase (left column) and the soft phase (right column) are shown separately. The value of , along the vertical axes, is normalized such that the area underneath the curve is unity. The results indicate that, throughout the deformation history, the strain is more strongly partitioned between the soft and the hard phase for the 2D microstructure than for the 3D microstructure (i.e. larger plastic strains, mainly in the soft phase). This explains why the macroscopic hardening mostly originates from the soft phase for the 2D microstructure. Furthermore, the plastic strain distribution is more heterogeneous for the 2D microstructure, because any local phase arrangement that promotes plastic straining extends infinitely in the outofplane direction (i.e. bands of soft phase at ). In the 3D case, the subsurface microstructure restricts the deformation, resulting in a less heterogeneous distribution of plastic strain.
A similar observation holds for the damage descriptor in the soft phase, shown in Figure 6(a–b). For the 2D microstructures, more grains have a higher value of . More importantly, also more grains display a value of , indicating fracture initiation. In other words, fracture initiation is ‘promoted’ in a twodimensional microstructure, for the same reasons as observed above.
4.3 Hotspot for fracture initiation
The fracture initiation hotspot is statistically identified by considering all unit cells. As discussed in Section 2.5, this leads to the probability of finding the hard phase at a certain position relative to fracture initiation: . The relative voxel coordinates are thereby aligned with the spatial coordinates . The result is shown in Figure 7 for the 2D microstructure and in Figure 8 for the 3D microstructure. Both of these figures include a crosssection along the plane. For the 3D microstructure a threshold is applied in order to visualize the most important morphological patterns. The probability scalebar is taken symmetric about the hard phase volume fraction to avoid an optical bias. The results may be interpreted as follows: (red) indicates an elevated probability of finding the hard phase at a position relative to the fracture initiation sites and (blue) an elevated probability of finding the soft phase there.
For the 2D microstructure, Figure 7 shows that, as assumed at the onset, fracture initiates in the soft phase which is indicated by the blue grain at . Directly next to fracture initiation in direction hard phase is found in most cases, since . Perpendicular, in direction, the soft phase is found with . Further away, in direction regions of the hard phase are found on both sides. The regions of the hard phase are interrupted by regions of the soft phase under degree angles with respect to the axis. This result is consistent with [12] even though a different constitutive model and damage indicator were used there.
For the 3D microstructure, the crosssection along the plane in Figure 8(a) closely resembles the result of the 2D microstructure. The probability is higher of finding hard phase () left and right, and soft phase () above and below. Figure 8(c) shows that fracture initiates where a ‘plate’ of hard phase is interrupted by two ‘channels’ of soft phase. Soft phase under degree angles is only recovered along the inplane direction.
5 Twodimensional vs. threedimensional deformation
The two different deformation paths from Section 2.3 are compared next, both applied to 3D microstructures.
5.1 Macroscopic response
The macroscopic responses are shown in Figure 9. The black curves correspond to applied planar shear and the green to applied axisymmetric shear, in between the responses of the individual phases (blue for soft and red for hard). There is little scatter, below 1%, for the different microstructures in the ensemble, i.e. the upper and lower bounds of all unit cells almost coincide. At the onset of yielding, slightly more hardening is observed under axisymmetric shear than for planar shear, as the green curves lie just above the black curves for (the difference is so small that it is difficult to see).
5.2 Microscopic plastic response
In addition to the response to planar shear, in Figure 4, the response to axisymmetric shear is shown in Figure 10. It is observed that the plastic strain response is more homogeneous for this case. Less extreme values are observed in the soft phase, and also the average of is smaller than for planar shear. At the same time the plastic strain is higher in the hard phase, where the average value is . In each of the three crosssections, small regions of elevated occur in the soft phase, at degree angles with respect to the three axes. These localized regions are located where soft phase is surrounded by hard phase. In addition, isolated grains are observed in the soft phase where is high. These regions are part of a previously described shear band along the inplane direction.
The probability density, , of the accumulated plastic strain, is shown in Figures 5(e–f) for axisymmetric shear. In comparing to the earlier results for planar shear (cf. to Figures 5(c–d)) it is observed that, throughout the deformation history, the plastic strain is more heterogeneous for planar shear than for axisymmetric shear. For the soft phase, Figure 5(d,f), a larger number of grains display a high plastic strain, and also the maximum is higher. For the hard phase, Figure 5(c,e), it is observed that the average plastic strain is higher for axisymmetric shear than for planar shear. Furthermore, the distribution is more skewed towards low values of for planar shear, indicating that most of the hard grains reveal a low .
5.3 Ductile fracture initiation
Figure 11 shows the damage descriptor, , in the microstructure of Figure 1 for an applied equivalent strain of for both load cases. The top row, Figure 11(a–c), shows three different crosssections parallel to the plane for the applied planar shear case. The bottom row, Figure 11(d–f), shows the same crosssections for the axisymmetric shear case. The different crosssections are taken at respectively one and two grains in negative direction with respect to the reference view on the left. The hard grains are colored gray, as damage in the hard phase is not considered.
For axisymmetric shear it is observed that the damage, , is lower in many of the grains when compared with the planar shear case (cf. Figures 11(a,d) for a single unit cell and Figures 6(b,c) for the ensemble). This is a direct consequence of the less heterogeneous distribution of the plastic strain. The purpose of the damage descriptor is however to identify fracture initiation, i.e. the focus should be on the sites where . The number of these fracture initiation sites is approximately the same (cf. Figures 6(b,c)). A large part of the fracture initiation sites coincide for the two load cases (e.g. the grain highlighted with a white circle), but also clear differences exist. One example is highlighted with a white square, for which the damage is high for axisymmetric shear, but not for planar shear (cf. Figures 11(a,d)). For this case, the morphology consisting of a soft grain flanked by hard grains on both sides in direction and soft grains on both sides in direction is not extended into the subsurface (Figure 11(e)). This shows that spatial phase distribution, directly adjacent to but also further away from the fracture initiation site, in combination with the deformation, has to be critical for fracture to initiate.
5.4 Hotspot for fracture initiation
The fracture initiation hotspot for axisymmetric shear is shown in Figure 12, with the same views as used in Figure 8. In the case of axisymmetric shear, the crosssection in Figure 12(a) reveals a pattern that is qualitatively the same as above. Adjacent to the soft grain in the center, in direction, hard phase is identified (), and in direction soft phase is found in most of the cases (). The influence of the morphology next to the fracture initiation site is thereby more pronounced than for planar shear. Further away, regions of hard phase are found in direction, which are interrupted by regions of the soft phase at degrees with respect to the axis. The latter are however less prominent than for planar shear. Extending the observations to three dimensions in Figure 12(b–c), the regions of hard phase in direction tend to be axisymmetric around the axis. Regions of the soft phase are found in cones at around the axis. Note that for all neighboring voxels at the resulting indicates soft phase, but for some it is below the threshold value.
In earlier work, De Geus et al. [12] used a simple mechanical analysis to confirm that the combination of a band of hard phase in the direction of applied tension, in combination with a phase boundary perpendicular to it, results in a hydrostatic tensile stress state. The regions of the soft phase in the directions of principal shear induce a large plastic strain. The combination of these morphological features promotes fracture initiation. This also explains the results for the threedimensional microstructure and deformation. For planar shear, in Figure 8, the plate of hard phase corresponds to a phase boundary whose area perpendicular to the tensile axis is maximized in such a way that the soft regions at degrees with respect to the axis are uninterrupted. The principal shear acts only in the plane and therefore soft phase is only recovered in that plane. For axisymmetric shear, the principal shear acts on all planes at degree angles around fracture initiation, explaining the resulting morphology of soft phase. The observed pattern of hard phase again maximizes the area of the phase boundary perpendicular to the axis of applied tension avoiding interruption of the shear bands.
6 Concluding remarks
A twophase microstructure is studied for which the two phases have different mechanical properties (yield stresses), resulting in a sharp mechanical contrast. Consequently, fracture initiation in the soft phase is strongly influenced by the spatial distribution of the phases. A large ensemble of random threedimensional microstructures is analyzed to systematically characterize the configuration that is most critical for fracture initiation. To this end, all fracture initiation sites are considered at once by calculating the average phase distribution around fracture initiation sites.
Fracture preferentially initiates when the following conditions apply:

regions of soft phase are aligned with the shear directions;

phase boundaries perpendicular to the direction of principal strain, i.e. hard phase on both sides, aligned with the tensile axis;

regions of hard phase are located next to the shear bands.
Using a twodimensional microstructure, the same conclusions are obtained. Indeed, the average phase distribution is a crosssection of the result obtained in 3D. However, the strain is more strongly partitioned between the phases as the deformation is not restricted or promoted by the subsurface microstructure. As a result, fracture initiation is more frequently observed in 2D at the same overall strain.
Acknowledgments
This research was carried out under project number M22.2.11424 in the framework of the research program of the Materials innovation institute M2i (www.m2i.nl).
References
plus 0.3ex
References
 Koo et al. [1980] J.Y. Koo, M.J. Young, and G. Thomas. On the law of mixtures in dualphase steels. Metall. Trans. A, 11(5):852–854, 1980. doi: 10.1007/BF02661217.
 Davies [1978] R.G. Davies. The deformation behavior of a vanadiumstrengthened dual phase steel. Metall. Trans. A, 9(1):41–52, 1978. doi: 10.1007/BF02647169.
 Maire et al. [1997] E. Maire, D.S. Wilkinson, J.D. Embury, and R. Fougères. Role of damage on the flow and fracture of particulate reinforced alloys and metal matrix composites. Acta Mater., 45(12):5261–5274, 1997. doi: 10.1016/S13596454(97)00147X.
 González and LLorca [2000] C. González and J. LLorca. A selfconsistent approach to the elastoplastic behaviour of twophase materials including damage. J. Mech. Phys. Solids, 48(4):675–692, 2000. doi: 10.1016/S00225096(99)000575.
 Brockenbrough and Suresh [1990] J.R. Brockenbrough and S. Suresh. Plastic deformation of continuous fiberreinforced metalmatrix composites: Effects of fiber shape and distribution. Scr. Metall. Mater., 24(2):325–330, 1990. doi: 10.1016/0956716X(90)90264H.
 Brockenbrough and Zok [1995] J.R. Brockenbrough and F.W. Zok. On the role of particle cracking in flow and fracture of metal matrix composites. Acta Metall. Mater., 43(1):11–20, 1995. doi: 10.1016/09567151(95)902562.
 Povirk [1995] G.L. Povirk. Incorporation of microstructural information into models of twophase materials. Acta Metall. Mater., 43(8):3199–3206, 1995. doi: 10.1016/09567151(94)004873.
 Schröder et al. [2010] J. Schröder, D. Balzani, and D. Brands. Approximation of random microstructures by periodic statistically similar representative volume elements based on linealpath functions. Arch. Appl. Mech., 81(7):975–997, 2010. doi: 10.1007/s0041901004623.
 Kumar et al. [2006] H. Kumar, C.L. Briant, and W.A. Curtin. Using microstructure reconstruction to model mechanical behavior in complex microstructures. Mech. Mater., 38(810):818–832, 2006. doi: 10.1016/j.mechmat.2005.06.030.
 LLorca and Segurado [2004] J. LLorca and J. Segurado. Threedimensional multiparticle cell simulations of deformation and damage in spherereinforced composites. Mater. Sci. Eng. A, 365(12):267–274, 2004. doi: 10.1016/j.msea.2003.09.035.
 Segurado and LLorca [2006] J. Segurado and J. LLorca. Computational micromechanics of composites: The effect of particle spatial distribution. Mech. Mater., 38(810):873–883, 2006. doi: 10.1016/j.mechmat.2005.06.026.
 de Geus et al. [2015] T.W.J. de Geus, R.H.J. Peerlings, and M.G.D. Geers. Microstructural topology effects on the onset of ductile failure in multiphase materials  A systematic computational approach. Int. J. Solids Struct., 6768:326–339, 2015. doi: 10.1016/j.ijsolstr.2015.04.035. arXiv: 1604.02858.
 Lewandowski et al. [1989] J.J. Lewandowski, C. Liu, and W.H.J. Hunt. Effects of matrix microstructure and particle distribution on fracture of an aluminum metal matrix composite. Mater. Sci. Eng. A, 107:241–255, 1989. doi: 10.1016/09215093(89)903924.
 de Geus et al. [2016] T.W.J. de Geus, R.H.J. Peerlings, and M.G.D. Geers. Competing damage mechanisms in a twophase microstructure: how microstructure and loading conditions determine the onset of fracture. Int. J. Solids Struct., Accepted, 2016. doi: 10.1016/j.ijsolstr.2016.03.029. arXiv: 1603.05841.
 Ramazani et al. [2013] A. Ramazani, P.T. Pinard, S. Richter, A. Schwedt, and U. Prahl. Characterisation of microstructure and modelling of flow behaviour of bainiteaided dualphase steel. Comput. Mater. Sci., 2013. doi: 10.1016/j.commatsci.2013.05.017.
 Xu et al. [2013] W. Xu, X. Sun, D. Li, S. Ryu, and M.A. Khaleel. Mechanismbased representative volume elements (RVEs) for predicting property degradations in multiphase materials. Comput. Mater. Sci., 68:152–159, 2013. doi: 10.1016/j.commatsci.2012.10.026.
 Moorthy and Ghosh [1998] S. Moorthy and S. Ghosh. A Voronoi Cell finite element model for particle cracking in elasticplastic composite materials. Comput. Methods Appl. Mech. Eng., 151(34):377–400, 1998. doi: 10.1016/S00457825(97)001606.
 Ghosh and Moorthy [2004] S. Ghosh and S. Moorthy. Three dimensional Voronoi cell finite element model for microstructures with ellipsoidal heterogeneties. Comput. Mech., 34(6):510–531, 2004. doi: 10.1007/s0046600405985.
 Moulinec and Suquet [1998] H. Moulinec and P. Suquet. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Methods Appl. Mech. Eng., 157(12):69–94, 1998. doi: 10.1016/S00457825(97)002181.
 Johnson and Cook [1985] G.R. Johnson and W.H. Cook. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Eng. Fract. Mech., 21(1):31–48, 1985. doi: 10.1016/00137944(85)900529.
 Sun et al. [2009] X. Sun, K.S. Choi, W.N. Liu, and M.A. Khaleel. Predicting failure modes and ductility of dual phase steels using plastic strain localization. Int. J. Plast., 25(10):1888–1909, 2009. doi: 10.1016/j.ijplas.2008.12.012.
 AlAbbasi and Nemes [2003] F.M. AlAbbasi and J.A. Nemes. Micromechanical modeling of dual phase steels. Int. J. Mech. Sci., 45(9):1449–1465, 2003. doi: 10.1016/j.ijmecsci.2003.10.007.
 Vajragupta et al. [2012] N. Vajragupta, V. Uthaisangsuk, B. Schmaling, S. Münstermann, A. Hartmaier, and W. Bleck. A micromechanical damage simulation of dual phase steels using XFEM. Comput. Mater. Sci., 54:271–279, 2012. doi: 10.1016/j.commatsci.2011.10.035.
 Brisard and Dormieux [2010] S. Brisard and L. Dormieux. FFTbased methods for the mechanics of composites: A general variational framework. Comput. Mater. Sci., 49(3):663–671, 2010. doi: 10.1016/j.commatsci.2010.06.009.