## 1 Introduction  
With the rapid development of modern industrial technology, particle-reinforced metal matrix composites, which combine the ductility and toughness of metal phase with the high strength of the reinforced particles, have been widely used due to their high specific strength and high specific rigidity [1; 2; 3; 4]. As is well known, the mechanical behaviors of composites with complex microstructures are closely dependent on many factors, such as inclusion volume fraction, inclusion shapes, etc. The higher stiffness behaviors of composites can be achieved by increasing the particle volume fraction (PVF) and improving the manufacturing technology [5; 6]. However, increasing the volume fraction often decreases the ductility of composites. Therefore, it is important for designers to consider some certain factors including the PVF, interphase between particles and matrix, and particle size when determining the mechanical properties. For estimating the particle size [7; 8; 9; 10] and particle shape [11; 12] effect on the mechanical properties, the representative volume element (RVE) is always used to establish the mechanical model [13; 14; 15; 16]. At present, a number of modeling techniques have been proposed. The strain-gradient plasticity models have been used by [PERSON] and [PERSON] [17; 18] to study particle size effect and interface fracture in aluminum alloy composites. [PERSON] and [PERSON] [19] modified the tangent-based homogenization method with incorporation of dislocation strengthening. It was revealed that the nonlinear behaviors of composites in uniaxial tension are based upon particle size and PVF. [PERSON] et al. [20] presentedan enhanced continuum model to explore the effects of particle size and volume fraction on the variability of the inelastic deformation response of particle-reinforced metal matrix composites. It was revealed that the variance of the composite response increases as particle size decreases and as the volume fraction of reinforcing particles increases. [PERSON] et al. [PERSON] et al. (2015) employed the finite element method to calculate the average coefficient of thermal expansion, [PERSON]'s modulus, and [PERSON]'s ratio. [PERSON] et al. (2016) combined the finite element method with statistical methods to study the elastic properties of composites with natural particle fillers. In order to investigate the effect of the interphase between the particles and matrix, a number of numerical techniques have been improved in the literature [PERSON] et al. (2015); [PERSON] et al. (2016); [PERSON] et al. (2017); [PERSON] et al. (2018). By incorporating the dislocation punched zone model, cohesive zone model, and the Taylor-based nonlocal theory, an enhanced finite element method (FEM) was proposed by [PERSON] et al. [PERSON] et al. (2016) to investigate the effect of interfacial debonding in the particle-reinforced metal matrix composites. [PERSON] et al. [PERSON] et al. (2018) combined the cohesive element technique with the finite element method to study the effect of the interface thickness and strength. It was shown that the particle shape and interface geometry strongly influenced the distribution of the stress states. Based on the coupled stress and energy criterion of the finite fracture mechanics (FFM), [PERSON] (2018) proposed a theoretical model to study the influence of the interface properties. From the investigations mentioned above, the researchers have always paid attention to the analysis of interfacial debonding effects on nonlinear stress-strain behaviors. [PERSON] and [PERSON] [PERSON] and [PERSON] (2018) employed the micromechanical FEM to predict the interphase damage and residual stresses effect on the mechanical behaviors of particle-reinforced composites (PRC). [PERSON] et al. [PERSON] et al. (2018) employed the commercial finite-element code ANSYS to establish a three-phase model consisting of a matrix, particles, and an interphase, in order to investigate the fracture behaviors of particulate polymer composites.