Statistical Analyses of the Strengths of Particulate Reinforced Metal Matrix Composites (PRMMCs) Subjected to Multiple Tensile and Shear Stresses

For design and application of particulate reinforced metal matrix composites (PRMMCs), it is essential to predict the material strengths and understand how do they relate to constituents and microstructural features. To this end, a computational approach consists of the direct methods, homogenization, and statistical analyses is introduced in our previous studies. Since failure of PRMMC materials are often caused by time-varied combinations of tensile and shear stresses, the established approach is extended in the present work to take into account of these situations. In this paper, ultimate strengths and endurance limits of an exemplary PRMMC material, WC-Co, are predicted under three independently varied tensile and shear stresses. In order to cover the entire load space with least amount of weight factors, a new method for generating optimally distributed weight factors in an n dimensional space is formulated. Employing weight factors determined by this algorithm, direct method calculations were performed on many statistically equivalent representative volume elements (SERVE) samples. Through analyzing statistical characteristics associated with results the study suggests a simplified approach to estimate the material strength under superposed stresses without solving the difficult high dimensional shakedown problem.


Introduction
Owing to its outstanding strength-to-weight ratio and excellent global isotropy, particulate reinforced metal matrix composites (PRMMCs) consisted of ductile metal matrix and discrete ceramic reinforcement particles became increasingly popular in recent years. In numerous applications where the weight of a structural component is a critical issue, besides optimizing the geometry and topology of a part, another plausible strategy is to replace the metallic materials by PRMMCs. To fully exploit the potential of PRMMC materials, it requires their mechanical behaviors as well as the relationships between these behaviors and underlying material structure to be well understood. Since experimental investigations alone provide only limited information on these issues, advanced micromechanical models were developed and employed as supplementary means. Using these models it was carefully examined and clarified in a large body of works how multiple microstructural characteristics, e.g. particle shape [1,2], size [3,4], distribution [5,6] individually and jointly influence global material properties ranging from coefficient of thermal expansion [7] to the rate of crack propagation [8].
Among global material behaviors investigated, the load bearing capacity of PRMMCs occupies a central place. While due to the complexity of service condition that most mechanical parts made by PRMMCs are subjected to, it is essential that the material strengths under both monotonic and cyclic loading can be accurately predicted. To this end, direct method based on Melan and Koiter's theorem has a remarkable advantage over the conventional incremental method since they allow the strengths to be calculated without considering the specific load history. Using direct methods (DM) constructed from static theorem in conjunction with the finite element (FE) based homogenization techniques, Weichert et al. [9], Magoariec et al. [10], Bourgeois et al. [11], Zhang et al. [12], You et al. [13], Chen et al. [14], Zhang et al. [15,16] elucidated in their works how the macroscopic feasible load domains of different heterogeneous materials can be calculated. Similarly, macroscopic strengths were also evaluated from numerical methods developed based on kinematic theorem, Carvelli [17], Chen and Ponter [18], Barrera et al. [19], Li and Yu [20,21], Canh et al. [22].
Strengths of PRMMC materials are not only affected by the material microstructure, but also by the type of loading they are exerted to. Compared to the uniaxial loading, microstructure plays a more important role in determining the global material strength under multiaxial loading featured by superimposed tensile and shear stresses [23,24]. When the material strength is predicted in such circumstance, i.e., given that multiple tensile and shear stresses vary independently, then in contrast to the severe difficulty by using incremental method the problem can be tackled more easily and with a more affordable cost by adopting direct method thanks to their path-independent nature. Extending direct method to high dimensional load space and improving the numerical efficiency in solving the resulting optimization problem have been extensively studied in several works, e.g., Kim and You [25], Simon and Weichert [26] explained how to formulate and solve such a problem by nonlinear programming, Spiliopoulos and Panagiotou [27] by a novel numerical procedure called Residual Stress Decomposition Method for Shakedown (RSDM-S), and Peng et al. [28] by a stress compensation method.
One intrinsic difficulty related to micromechanical analysis of PRMMC materials is that, unlike most fiber reinforced composites, the material often has a random microstructure. As pointed out in many research works [7,29,30], predicting material behaviors of PRMMCs from unit cell models or representative volume elements (RVE) having an improper size may lead to strongly biased results. To account for the randomness in the material microstructure and to objectively reflect its influence in the model, one generally acknowledged solution is to employ statistically equivalent representative volume elements (SERVEs) [31]. According to this theory, instead of from one individual RVE, the material behavior should be viewed from a statistical point of view and predicted from a set of SERVE models. One key premise for studying PRMMC materials by SERVEs is to have capable algorithms for creating a sufficient number of samples. In this regards, algorithms such as random sequential absorption and modeling software such as Digimat and Dream3D have greatly reduced the effort required for generating samples and therefore make this strategy a viable means for studying PRMMC materials.
In our previous works [32,33], an integrated numerical approach for predicting the ultimate strength and endurance limit of PRMMC materials is proposed which incorporates homogenization, direct methods, and statistical analysis. In the present study, this method is extended to the case of high dimensional load space such that the strengths of PRMMCs subjected to superimposed tensile and shear stress can be evaluated. To exhaustively investigate the material behavior under all plausible combinations between two tensile and one shear stresses, a method for generating optimally distributed weight factors in an n dimensional space is developed. Using this method, endurance limit of an exemplary PRMMC material, WC-30 Wt.% Co, submitted to three stresses is evaluated from many SERVE samples.

Overall Strengths of Heterogeneous Materials
Macroscopic strengths of heterogeneous materials can be predicted by applying direct method to RVE samples. To do this, the material should be considered in two well-separated scales: The microscopic scale y in which the RVE exists and the macroscopic scale x large enough that the heterogeneities smear out. Two scales are related by a scale parameter δ For a heterogeneous material, once it is subjected to external loading, microscopic stress field σ and the macroscopic one Σ satisfy the relationship (2) Σ = 1 Ω � σ (y)dV = σ (y) . Here �·� stands for the mean-field averaging operator, and Ω indicates the RVE domain. Similarly, the relationship between strains in micro and macro scales satisfy The microscopic strain ε can be decomposed into two parts: The average value E and a fluctuating part ε * When the overall behavior of the RVE is purely elastic, then Σ and E are correlated by an effective elastic tensor C If the material is isotropic in the macro scale, then C can be uniquely determined from effective Young's modulus Ē and effective Poission's ratio ν.
Given that the heterogeneous material is composed of elasto-plastic constituents, its global ultimate strength U and endurance limit ∞ correspond to plastic and shakedown limits, respectively. Here, plastic limit can be viewed as a particular case of shakedown in which the load is only allowed to vary monotonically. The feasible load domain of a material can be calculated from DM by imposing either E or as boundary conditions [9], where the former case refers to the kinematically uniform boundary condition (KUBC) and the latter case the statically uniform boundary condition (SUBC).
Here, ∂� refers to the RVE surface, n the outer surface normal, and u * the fluctuation part of the displacement corresponds to ε * . Because the material studied in this work is non-periodic, some modifications were made on conditions given in Eqs. (6) and (7): Instead of enforcing the node-wise anti-periodicity of the stresses and periodicity of the fluctuating strain, it is only required �ε * � and �ρ� to be zero. In other words, when an RVE is subjected to SUBC, �ρ� is only required to be divergence free inside and have no contribution on the macroscopic stress.

Generating Optimally-Distributed Weight Factors in an n Dimensional Load Space
In the present study, the aim is to evaluate strengths of RVEs simultaneously subjected to two tensile and one shear stresses, all vary independently. To this end, we introduce in this section a pragmatic algorithm which systematically generates optimally-distributed weight factors in an n dimensional load space. When an elastic-perfectly plastic structure is submitted to NL independently varying loads P n (x, t) , P n can be decomposed into two parts: A time independent base P 0n (x) corresponds to the load pattern and an associated coefficient µ i describes the load magnitude. This way, as shown in Figure 1, a specific load history H can be defined as where H is a trajectory in the n dimensional load space L For a periodic heterogeneous material in the absence of the body force subjected to SUBC, Magoariec et al. [10] proposes a static shakedown condition which requires the stress field associated with the reference elastic body B E , σ e , and the time invariant residual stress field ρ to satisfy (6) σ e : anti-periodic on ∂� , u * periodic on ∂� , �σ � = � .
anti-periodic on ∂� . �ρ� = 0 . Here, µ − n and µ + n correspond to the lower and upper bound of a loading P n , respectively. L is a convex domain with NV = 2 n vertices. We denote each vertex by P k , so k could be any integer between 1 and NV.
Due to the linearity of elastic stresses, σ e can be seen as a superposition As Eq. (10) indicates, for a specific load history, only the magnitudes of loads are varying, thus it is justified to completely focus on coefficients µ n . By separating µ n from their original basis P 0n and merging with the unit basis e n , a load magnitude vector µ can be formed. All It is evident that an one-to-one mapping can be established between elements in L and U . The benefit for introducing U is, that the physical nature of a specific loading becomes unimportant; loadings P n differ only by their magnitudes. For simplicity, each basic loading P 0n should be adjusted to a comparable level, e.g. all correspond to the elastic limit. Then with µ n varying inside µ n ∈ [µ − n , µ + n ] , all admissible load combinations can be exhausted.
To illustrate the algorithm for generating µ with the lowest complication, the discussion is restricted to a particular load case where any given load P n is forced to vary within [0, µ + nP 0n ] . In spite of this restriction, it is clear that by replacing 0 in the n th entry in result by µ − nP 0n , the outlined approach can be easily generalized to account for reversed loading as well. To derive σ e at all NV load vertices in an n dimensional load space, we define a matrix U ∈ R NV ×n whose each row corresponds to a binary index of a vertex k Based on U , a load magnitude vector µ can be converted to a combinatorial matrix U µ ∈ R NV ×n (10) The purpose of U µ is to superpose stress calculated from each individual loading P 0n to all NV vertex stresses σ e k . To do this, first we arrange stress associated with each basic loading P 0i in following matrix form In this equation, σ e (P 0n ) is a column vector consisting of stresses at all NG Gaussian points. When a 3D model is considered, the length of σ e i (P 0n ) is 6 · NG . Using σ e B introduced in Eq. (14), σ e k at all load vertices can be calculated by In DM, load factor α is often evaluated for many different combinations of P i . As Eq. (13) implies, each combination can be viewed as a polyhedron whose vertices generated from a different µ . A challenge in DM is to use the least number of µ to approximate the domain U . For two independently varied loads P 01 and P 02 , one simple way to generate a series of µ is introducing an angle θ ∈ [0, 2π] : (16) µ = cos θ sin θ T . Similarly, when three independently varying loads P 01 , P 02 and P 03 are considered, one can use two angles θ ∈ [0, π] and ϕ ∈ [0, 2π] Although this approach is simple, it is hard to be generalized to cases n > 3 . Meanwhile, the method requires users to specify the values of parameters, e.g., θ and φ in Eq. (17), therefore when these parameters are inappropriately selected, µ 's will have a severe inhomogeneous distribution in U . To remedy such problem, we developed a simplealgorithm. Before illustrating this algorithm, we assume that P 0n has already been adjusted to a comparable level, in other words, σ e (P 0n ) are in similar magnitudes. Based on this premise, we could fix all µ + n to 1 and discretize the interval [0, 1] associated to P 0n to a grid vector p n of (17) µ = sin θ cos ϕ sin θ sin ϕ cos θ T .
the length l n . Because each p n can be viewed as a set {p n } whose elements are vector's components, the load set U can be generated by calculating the Cartesian power p 1 × p 2 × · · · × p n . The U generated in such a manner is a set consisting of n i=1 l i elements where each element represents a vector µ ∈ R n . It is easy to figure out that many µ in U have different norms but same directions. These vectors are repetitive because their difference in norms can be compensated by the load factor α calculated from DM. Due to this reason, only vectors differing in terms of directions should be remained in U and the repeating ones should be removed. To do this, as shown in Figure 2 we normalize vectors in U and remove the repeating ones after normalization. It is noteworthy that, after doing this µ in U are no longer evenly distributed. To optimize the distribution, the distance between each paired vector µ is evaluated and once µ drops below a user-defined tolerance ǫ , two vectors are regarded as identical and one of them is removed from the set (see Figure 2(b)). Using this method, by adjusting parameter ǫ , one can easily control the size of U and generate the set with optimally distributed µ's.

Optimization Problem from the Static Shakedown Theorem
Employing shakedown conditions outlined in Eqs. (6) and (7) and discretizing the physical fields by FE formulations, the macroscopic strengths of RVEs composed of elastic perfectly plastic materials can be calculated from following optimization problem Here, α is the load factor, C the equilibrium matrix, ρ i the stress tensor in the i th Gaussian point, σ e ik the abbreviation of σ e i (P k ) which means σ e at Gaussian point i and load vertex k, σ Y the yield strength, F the yield function. Solving Eq. (18) yields the load capacity of the RVE, and depending on if k = 1 or k > 1 the calculated strength corresponds to either ultimate strength U or endurance limit ∞ .
Following steps elucidated in our previous work [32] problem Eq. (18) can be transformed into an equivalent one whose inequality constraints are Euclidean ball constraints. Since after this transformation the problems obtained are typical second-order cone programming (SOCP) problems, they can be solved efficiently using commercial optimization solvers such as Gurobi [34], CPLEX [35], MOSEK [36], among others. In the present study, solver Gurobi is employed to calculate strengths (18)  (b) After removing repeated ones

Configuration of FE Models
In the present study, an exemplary PRMMC material, WC -30 Wt.% Co, was chosen, and by adopting the numerical approach explained above we studied its mechanical behaviors under two tensile and one shear stresses. Microstructures of WC -30 Wt.% Co are similar to the scanning electron microscope (SEM) image shown in Figure 3 where the dark grey areas are Co and the bright areas are the WC-grains with their characteristic prismatic shape. The average size of carbide phase is d WC = 2.35 µ m. Using a modeling technique developed in-house, fifty SEM images were cropped to 40 µ m × 40 µ m squares and converted to FE models. All these models share a same mesh pattern: Elements distant from the phase boundaries were assigned a global size of 0.8 µ m and elements close to the phase boundaries with an edge length of 0.2 µ m. Meanwhile, as has been highlighted in our previous work [32], on the one hand, 3D models demand a huge amount of computational power to be calculated and therefore are too expensive for a statistical analysis, while on the other hand both plane strain and plane stress models make extreme idealizations of the material behavior in z direction and show a nontrivial mesh sensitivity. For these reasons, due to its limited modeling complexity and insignificant mesh size dependence, 2.5D models obtained by extracting 2D models for 1 µ m in the z direction were employed for the present study. Meanwhile, the study is restricted to the infinitesimal strain theory and consider both phases as von Mises materials with parameters summarized in Table 1.
To perform shakedown analysis and determine overall strengths of SERVE samples, elastic analyses were performed in the commercial FE code ABAQUS. In these analyses, as illustrated in Figure 4 all fifty samples were successively prescribed with global stresses 11 , 22 and τ 12 by means of SUBC and the elastic stress field under each load case was evaluated. Employing these stress fields and other key model information such as the shape function and node coordinates we constructed the shakedown problems Eq. (18) at 44 load combinations whose weight factors are calculated from the aforementioned algorithm. By solving these problems considering one, two or eight load vertices, for each sample three feasible load domains corresponding to U , 2P ∞ and 8P ∞ were obtained.

Statistical Analysis of Elastic Properties
Before analyzing the strengths of RVE samples, we first evaluated results gained from elastic analyses. The purpose for studying elastic behaviors is twofolds: Firstly, a careful and thorough investigation of elastic responses of SERVE samples to the global shear stress is the premise for understanding their performances, especially the load bearing capacities, under shear and superposed tensile and shear stresses. Secondly, in reality WC -30 Wt.% Co is an isotropic material, which means it has only two independent elastic parameters. Therefore, by comparing the effective shear modulus obtained from the mean-field homogenization Ḡ with the one calculated from Ē and ν and denoted as Ḡ * , one may know the degree of sample anisotropy and infer from this value if they satisfy basic model requirements. Figure 5 shows the statistical distribution of carbide volume fraction WC Vol.%, Ē and Ḡ . In this figure r stands for the Pearson correlation coefficient and the proximity of r to one implies that both Ē and Ḡ are determined by the carbide content. Result of the comparative study between shear moduli obtained from two distinct approaches is illustrated in Figure 6. In this figure |δ| is the relative discrepancy defined as |δ| = |Ḡ * −Ḡ|/Ḡ . Figure 6 shows that the difference between Ḡ * and Ḡ is insignificant for all SERVE samples. The fact that |δ| is no greater than 2.0% for all samples implies that all these models demonstrate a sound global isotropy. For this reason, it is justified to conclude that they satisfy the basic size requirement and therefore can be used to study more sophisticated nonlinear material behaviors.

Statistical Analysis of Composite Strengths
The average feasible load domain calculated from fifty SERVE samples can be found in Figure 7. In this figure all macroscopic strengths were normalized with respect to the binder yield limit in order to emphasize the strengthening effect of carbide particles. One can easily notice from this figure that in both � 11 − τ 12 and � 22 − τ 12 planes, 2P ∞ overlaps with 8P ∞ . The reason behind this phenomenon can be explained as follows: Failure of all samples is caused by the stress state at P 1 in Figure 1, since this vertex plays a dominant role, as a consequence it does not make any difference if the tensile and shear stresses have to be proportional or independently. The slight difference between 2P ∞ and 8P ∞ is attributed to the fixed solver configuration parameters: The optimization problem Eq. (18) considering two and eight load vertices are remarkably different in terms of both variables and constraints, hence it is unavoidable that solutions for two cases are slightly different when adopting the same termination criterion. In the 11 − 22 plane 2P ∞ is remarkably smaller than 8P ∞ , especially in the region close to the π/4 line. As explained in Ref. [32], this phenomenon is caused by the fact that when 11 and 22 vary proportional, they result in a hydrostatic stress which does not directly contribute to the overall strength.
Besides the mean value, the scatter of strengths among SERVE samples were evaluated as well and summarized in Figure 8. In this figure the interval between [−2SD, 2SD] where SD stands for the standard deviation is illustrated as a colored band. Since width of a colored band reveals how great the impact of a microstructure is, it is clear that with the increase of NV the influence of microstructure becomes less critical. The relationship between volume fraction and strengths of different kinds can be found in Figure 9. One can notice from this figure that, unlike elastic modulus, strengths of the PRMMC samples are not strongly correlated to the carbide volume fraction. Meanwhile, one can see that r between U in 11 and 22 directions is 0.47, while r between ∞ in two direction becomes 0.70. It means, compared to ultimate strengths, endurance limits of samples subjected to normal stresses in 11 and 22 directions are apparently more similar. This again confirms the finding that the microstructure has a greater influence on the material strength in the context of a monotonic load than a time-varied load. Another noteworthy phenomenon in Figure 9 is that for both U and ∞ , the correlation between tensile and shear stresses is insignificant. This implies that a microstructure results in an outstanding tensile strength does not necessarily results in an outstanding shear strength. To understand what are the main factors determining the shear strength of the PRMMC material, we compared stress fields obtained from the DM calculation between SERVE samples having poor and excellent shear strengths. As can be noticed from Figure 10, the distribution pattern of total stress σ = ασ e +ρ at vertex P 1 are remarkably different for two SERVE samples. For the sample having a poor shear strength whose Z-score, Z τ ∞ , is merely −2.82, the distribution of normalized effect stress is quite localized. It means that the overall deformation concentrates in a weakest region consisted almost completely of the binder phase when the sample  is subjected to a global stress. Undoubtedly for a sample whose deformation is dominated by such a mechanism, the reinforcement particles can hardly contribute and, as a consequence, the overall strength of the sample is poor. In contrast to these samples, for models in which such a weakest region is not obvious, a higher fraction of the overall volume would be stressed. Under such circumstance, the carbide particles would come into effect and the sample would have a higher strength.
To further investigate how SERVE samples react to superposed tensile and shear stresses, we separated the loads to two groups and considered shear stress as an addition to two tensile stresses. Adopting this point of view, we viewed 3D feasible load domain by means of contour maps. As shown in Figure 11, each contour line in such a map corresponds to a 2D load domain with a shear stress whose magnitude is displayed. Figure 11(a) corresponds to Figure 7 viewed from above, and Figure 11(b) shows load domains for a randomly picked sample. Contour maps reveal how strength under two tensile stresses change in accordance to the additionally imposed shear stress. These figures can be interpreted as a trade-off between two tensile and one shear stresses. One noteworthy characteristic pertained to Figure 11 is that the transition of contour lines alongside increase of τ is fairly steady: The shape of contour lines are well preserved, and the rate of their shrinkage for all combinations of 11 and 22 are alike. Such a characteristic can help with simplifying the strengths prediction practice of materials similar to WC -30 Wt.%. One should bear in mind that although by adopting the numerical method elaborated in the present work the entire feasible load  domain under three combined tensile and shear stresses can be calculated, the required computational power for solving shakedown problems in 3D load space is still much greater compared to 2D. For this reason, when instead of the accurate value it is only important to know how the load domain varies with the additional shear stress, then the effort for solving DM problems with eight vertices can be saved and one can simply calculate three 2D load margins and extrapolate them to a 3D surface.

Conclusions
In this paper, the numerical approach for predicting the strengths of PRMMCs is extended to take into account of independently varied two tensile and one shear stresses. A general algorithm is proposed to generate optimally distributed weight factors in an n dimensional load space. Based on this algorithm the study calculated ultimate strengths and endurance limits of fifty 2.5D SERVE samples modeled from real SEM image of WC-30 Wt% Co. Through studying the statistical characteristics of results, the study shows that shear strength is only subtly related to both tensile strength and carbide contents, while the magnitude of the shear strength is determined by the highly stressed volume fraction. In addition, by investigating the load domains by means of contour maps the study suggests to decouple the superposed loads and simplify a shakedown problem in 3D load space to three problems in 2D load space. In the present study, the failure mode is restricted to the plastic failure and defects such as micro pores introduced during the sintering process are neglected, the influence of these factors on material strength will be investigated in our future studies.