Packing and size effects in elastic-plastic particulate composites: Micromechanical modelling and numerical verification

The issue of applicability of the Morphologically Representative Pattern (MRP) approach to elastic-plastic composites is addressed. The extension to the regime of non-linear material behaviour is performed by employing the concept of incremental linearization of the material response in two basic variants: tangent and secant. The obtained predictions are evaluated through comparison with the outcomes of numerical analyses. Finite Element simulations are carried out using periodic unit cells with cubic arrangements of spherical particles and representative volume elements (RVE) with 50 randomly placed inclusions. In addition to the analysis of the packing effect in two-phase composites, the size effect is also studied by assuming an interphase between the matrix and inclusions. It is concluded that the MRP approach can be used as an effective predictive alternative to computational homogenization, not only in the case of linear elasticity but also in the case of elastic-


Introduction
Particulate composites often exhibit packing and size effects. For example, the influence of particle packing can be observed when the presence of agglomerates deteriorates a composite's fracture toughness ( Kushvaha & Tippur, 2014 ) or when the yield stress increases with a decreasing distance between particles ( Kouzeli & Mortensen, 2002 ). As concerns the size of constituents, its basic impact is related to the role of interfaces between the phases, e.g. due to chemical reactions between metallic and ceramic materials or special treatment of a particle's surface, this zone can have different properties than the two basic components ( Chmielewski et al., 2017 ). Moreover, in this particular place, phenomena degrading a composite's properties usually start, cf. ( Jarz ąbek, 2018 ). One may indicate in this respect debonding, cracking or damage evolution. These phenomena are strictly connected with stress and strain distribution in the heterogeneous material during the deformation process. Therefore, accounting for distinct interface properties when predicting the overall response can enhance the modelling accuracy.
The main issue that needs to be addressed in both approaches is the specification of the material properties of a 2D interface or the interphase, which in the latter case can be uniform or vary along the thickness as in Shen and Li (2003) , Sevostianov and Kachanov (2006) .
The focus of the present study is on the morphologically representative pattern (MRP) approach. For elastic composites, extensive verification of the MRP approach was conducted by Majewski et al. (2017) , where MRP estimates were compared with results obtained by the finite element computational homogenization. The present contribution continues this research and addresses the issue of applicability of the MRP approach to elastic-plastic composites. As it is usually done in the literature related to mean-field micromechanical models, an extension is performed by employing the concept of incremental linearization of the material response ( Hill, 1965;Kursa, Kowalczyk-Gajewska, Lewandowski, & Petryk, 2018;Ponte Castañeda & Suquet, 1997;Suquet, 1995;Tandon & Weng, 1988 ). Two basic variants of such linearization are tried out: tangent ( Hill, 1965 ) and secant ( Tandon & Weng, 1988 ). The obtained predictions are evaluated through comparison with the outcomes of corresponding FE analyses. FE simulations are carried out using periodic unit cells with cubic arrangements of spherical particles and representative volume elements (RVE) with 50 randomly placed inclusions. In the latter case, the procedure for generating RVEs which allows more control over the resulting mean minimum distance between particles (a packing parameter), presented in ( Majewski et al., 2017 ), is applied. Besides the impact of packing on the response of two-phase composites, also the size effect is investigated by assuming an interphase between the matrix and inclusions.
The paper is constructed as follows. In the second section, the formulation of the MRP mean-field model of composite materials described by a linear constitutive law in the phases is shortly recalled and then its extension to the elastic-plastic non-linear case is proposed. Section 3 describes the details of the numerical finite element analyses and the procedures used for generating periodic unit cells and representative volume elements with random particle distributions of varying values of the packing parameter. In Section 4, the predictions of the composite response in the elastic-plastic regime obtained Fig. 1. The microstructure of a two-phase composite reinforced with spherical inclusions k = 1 , . . . , 9 of the same radius R i is represented in the MRP model by two patterns: 4-GSC and SC. Above, λ k is the distance between inclusion k and its nearest neighbour, λ is the average value of λ k for all inclusions in the representative volume, which is λ = 1 9 9 k =1 λ k in the provided example, t int is the interphase thickness, and HEM is the Homogeneous Equivalent Medium.
by FE analysis and micromechanical modelling are compared. Finally, the last section includes a summary and the main conclusions.

Linear constitutive law
The Morphologically Representative Pattern (MRP) approach ( Bornert et al., 1996;Marcadon et al., 2007 ) is based on the micromechanical model of the composite sphere formulated by Hashin (1962) for a two-phase material, and its further modification towards the 3-phase Generalized Self-Consistent (GSC) scheme by Christensen and Lo (1979) , Christensen (1990) . In the following years, the GSC model was extended from the three-phase to the n -phase configuration by Herve and Zaoui (1993) , who solved the problem of the composite sphere with n -coatings and formulated the n -phase Generalized Self-Consistent (n-GSC) scheme.
The MRP model, set within the mean-field modelling framework, postulates a subdivision of the Representative Volume Element (RVE) of volume V with some morphological characteristics into M representative patterns α. In the frame of the model, under the assumption of a linear constitutive law, the linear relationship between the auxiliary far-field strain E 0 and the average strain ε α k in the phase k and the pattern α is described by the fourth order concentration tensor A α k , that is Then, observing that the overall average strain of the composite E equals ε V , where · V is defined as the volume averaging and c α is the volume fraction of the pattern α in the representative volume V , while f α k is the volume fraction of the phase k in the subvolume V α occupied by the pattern α. Utilizing the constitutive law σ α k = L α k · ε α k and the expression for the overall average stress in which the fourth order tensor L is the overall composite stiffness.
In this paper, the basic configuration of the MRP model is explored. The composite microstructure is represented using two simple patterns ( M = 2 ): 3-or 4-GSC and SC ( Fig. 1 ) which, despite the simplicity of the configuration, allow one to study the influence of size and packing effects. The morphology of the first pattern, the 4-GSC scheme, is presented in Fig. 1 b: it contains a composite sphere made of a spherical inclusion of radius R i , an inner coating of the interphase of thickness t int , and an outer coating of the matrix phase of thickness λ / 2 . The composite sphere is embedded in a Homogeneous Eqiuvalent Medium (HEM) with effective properties. The second pattern has a self-consistent type morphology with an inclusion made of pure matrix material ( Fig. 1 c). The thickness t int in the first pattern is independent of the particle radius R i and the interphase has different properties than the basic two phases: matrix and inclusions, so that t int accounts for the size effect . The parameter λ , in turn, accounts for the packing effect . In such a 2-pattern approach, the volume fractions of the GSC-type pattern, c 1 , and the SC-type pattern, c 2 , are calculated as where f i is the volume fraction of inclusions in the composite, R i is the particle radius, and λ is the mean distance between the nearest-neighbour particles. It is assumed that the composite components have isotropic properties, that inclusions are spherical and that they are isotropically placed in the matrix. It is easy to conclude that on such assumptions the concentration tensors A α k are isotropic fourth-order tensors specified by two scalars: α α P k and α α D k which depend on the pattern morphology, namely: I and I are the identity tensors of the fourth and second rank, respectively. In the present study, there are two basic schemes taken into account: the morphology of the n -phase model ( Fig. 1 b), and the simple morphology of the Eshelby problem ( Fig. 1 c). The relevant formulae for the concentration factors α α P k and α α D k can be found, e.g., in ( Herve & Zaoui, 1990;1993;Majewski et al., 2017 ).
Let us recall that, because the MRP model combines the n -GSC and SC schemes, the MRP estimates of effective properties for the analysed case fill the space between the GSC and SC estimates, and for a specified f i the MRP values depend on the packing parameter λ / ( 2 R i ) . Two limit cases of the packing parameter λ / ( 2 R i ) can be distinguished: • all the matrix material is used to form the coating of the inclusion in the GSC pattern: • each inclusion touches its nearest neighbour: λ / ( 2 R i ) = 0 .
In the first case, the volume fraction of the SC pattern is zero, so the MRP solution reduces to the GSC estimate. In the second case, the GSC pattern of the MRP model turns into the classical SC model and the MRP estimate coincides with the SC solution.
In Majewski et al. (2017) , the MRP model of the present form was thoroughly discussed and verified against FE calculations for elastic composites. The limit cases of infinitely stiff phases and voids were also considered. In this work, the framework is extended to elastic-plastic metal matrix composites and subsequently verified with respect to FE calculations performed on unit cells with regular arrangements of particles and representative volume elements with random particle distributions.

Extension to the elastic-plastic response
In order to take advantage of the micro-macro mean-field models, outlined in the previous section in the context of linear properties, incremental linearization proposed by Hill (1965) is adopted for elastic-plastic materials. Depending on the definition of the current stiffness tensor, two variants of the linearization procedure: tangent ( Hill, 1965 ) and secant ( Tandon & Weng, 1988 ) are considered. At each strain increment, the linearized constitutive relations in the phases for the tangent (superscript t) and secant (superscript s) method are: respectively, where the current tangent ( L t ) or secant ( L s ) elastic-plastic stiffness tensor is applied. The actual specification of these fourth order tensors depends on the constitutive model of the phases. Similarly to Kursa, Kowalczyk-Gajewska, and Petryk (2014) , Kursa et al. (2018) , the effective elastic-plastic response of metal-ceramic composites (MMC) with an elastic-plastic matrix and elastic inclusions is analysed. The inclusion phase is assumed to be a linear elastic material with parameters relevant to ceramics ( Table 1 ). The metal matrix is assumed to be made of a ductile material governed by linear elasticity and the Huber-von Mises yield function f ( σ m ) with the associated flow rule: where λ is a non-negative plastic multiplier ( λ ≥ 0) and σ m is the deviatoric part of the stress tensor. Additionally, the plastic multiplier and the function f ( σ m ) fulfil the consistency conditions: The isotropic hardening in the form of a power law is assumed for the yield stress Y ε p eq : Assumed material parameters of the metal matrix composite with an elastic-plastic matrix and elastic inclusions ( Suquet, 1997  The material parameters of the Metal Matrix Composite (MMC) are listed in Table 1 .
Using equations ( 7) -(10) , the tangent elastic-plastic stiffness tensor L t m in Eq. (6) a can be expressed by the following formula: in which: K m , G m are the elastic bulk and shear modulus, respectively. It should be noted that, in spite of the assumed isotropy of the constituent phases, the current tangent stiffness matrix L t m is anisotropic ( Chaboche et al., 2005;Doghri & Ouaar, 2003 ). To avoid excessive stiffness of the elastic-plastic response of the composite ( Chaboche et al., 2005;Doghri & Ouaar, 2003 ), as well as to be able to use the formula for concentration tensors for the n -GSC model available only for isotropic materials ( Herve & Zaoui, 1993 ), the current elastic-plastic stiffness tensor of the metal matrix obtained with the tangent incremental variant is subjected to the isotropization according to Chaboche et al. (2005) , Sadowski, Kowalczyk-Gajewska, and Stupkiewicz (2017) , Kursa et al. (2018) , namely Following Kursa et al. (2018) , in the case of secant linearization of the constitutive equations, the current elastic-plastic stiffness tensor of the matrix, for a proportional loading path, is defined in the isotropic form: The non-linear elastic-plastic response of a composite is obtained on the basis of the mean-field MRP averaging scheme described in the previous section in a step-by-step manner. The current stiffness tensor, determined by Eqs. (13) and (14) for the tangent and secant approach, respectively, is updated with the current accumulated plastic strain ε p eq . In the tangent variant, the incremental stress and strain are used in the algorithm, while in the secant approach, the current stress and strain tensors are sufficient to derive L s m . The algorithm used in this work for the MRP model is based on the previously developed algorithm for elastic-plastic composites. The previous algorithm, which computes the overall response of an elastic-plastic composite using an incremental iterative procedure, has been presented by Kursa et al. (2018) in Box 1. For the purposes of the present work, a modification is introduced which allows one to calculate MRP concentration factors for the composite phases at each incremental time step t + t. Concentration factors are calculated for the elastic inclusion phase, the elastic-plastic interphase, and the matrix and HEM phases. The computer code has been implemented in the Wolfram Mathematica ( www.wolfram.com ) environment.
The basic features of the two linearization methods when combined with the MRP micromechanical model are demonstrated in Fig. 2 . The responses of the individual components of the metal matrix composite, listed in Table 1 , are marked in gray colour in both graphs in this figure. The estimated response of the composite in isochoric extension (cf. Eq. (21) a), depending on the adopted linearization variant, is presented in Fig. 2 a. The composite response obtained using secant linearization (subscript S) is stiffer than the one from tangent linearization (subscript T). For the same linearization scheme, the estimate obtained with the MRP model falls between the SC and GSC estimates, consistently with the previous results regarding elastic properties of composites ( Majewski et al., 2017 ). Fig. 2 b contains results obtained using tangent linearization combined with different mean-field models: the self-consistent (SC), the generalised self-consistent (GSC), the Mori-Tanaka (MT), and the morphologically representative pattern (MRP) approach. The MRP estimates depend on the packing parameter. When the packing parameter tends to zero, the MRP approach delivers the same results as the SC scheme. Conversely, when the packing parameter tends to the maximum value equal to 0.489 for f i = 30% , then the MRP estimates overlap with the predictions of the GSC scheme.

Numerical models
Procedure for generation of cubic volume elements with randomly-placed inclusions. To test the influence of sphere packing on the material response more thoroughly, simulations were performed on samples with both regular and random arrangements of inclusions. Simulation results for composites with regularly-placed inclusions are presented in Section 4.1 , while those for metal matrix composites with random arrangements of inclusions are described in Section 4.2 .
The random samples-cubic volume elements with randomly-placed inclusions-were generated using a discrete element method-based dynamic procedure. This procedure was chosen over other, mainly static, available placement procedures because of its ease of implementation in existing software, speed of execution, and ability to handle large volume fractions (up to around 0.65). The procedure itself has been more thoroughly described in the previous work ( Majewski et al., 2017 ). Below we only summarize its most important features.
The generated microstructures occupy a 1 × 1 × 1 periodic cube, and are produced with the aid of the discrete element method system Yade ( Šmilauer et al., 2010 ) using a technique similar to that employed by Zieli ński (2015) , as follows. First, n non-overlapping elastic frictionless spheres of radius R i are placed in an enlarged cube of dimensions 2 × 2 × 2. Together, they take up a prescribed volume fraction f i of the target unit cube, i.e. n 4 3 π R 3 i = f i . In fact, the spheres are artificially increased by the factor d = 0 . 001 prior to the placement, to ensure that there remains a gap between the actual inclusions of radius R i . This initial, sparse random placement of spheres is obtained using a built-in random sequential addition ( Torquato, 2002;Widom, 1966 ) functionality of Yade. After the preparatory stage, the periodic cell shrinks uniformly to its target dimensions 1 × 1 × 1, pushing the spheres closer together. During the compression, the spheres mix by undergoing elastic collisions and crossing the periodic walls. Once the cell reaches unit size, a check is performed. If the sphere centers are separated by more than (2 + d) R i , i.e., the enlarged spheres are compressed by no more than d R i so that there remains a gap of at least d R i between the corresponding spheres of radius R i , then the positions of the spheres are adopted as the positions of inclusions in the composite.
For any given values of n and f i , the basic procedure generates random microstructures whose mean nearest neighbour distance λ is not controlled-it is an unpredicted outcome of each run of the procedure. To have more control over λ , two modifications were introduced. The first one consists in artificial enlargement of the spheres' radius, which forces computed sphere centers to lie farther apart-thereby increasing λ . The second one consists in using pairs of spheres kept a fixed distance apart instead of single spheres-this ensures that λ is not greater than the gap which is the same between the spheres of each pair. These two methods allow easier generation of microstructures with low and high values of λ . Example microstructures differing in their values of λ are shown in Fig. 3 .
Numerical homogenization procedure. To investigate numerically the response of the material under an applied global deformation, periodic statistical volume elements ( Ostoja-Starzewski, 2006 ) containing 50 inclusions each were adopted as approximate representative volume elements (RVE). Smaller samples with 30 inclusions have been previously used, e.g., by Pierard, LLorca, Segurado, and Doghri (2007) , Czarnota, Kowalczyk-Gajewska, Salahouelhadj, Martiny, and Mercier (2015) , to study properties of elastic-viscoplastic composites. Such samples are generally considered to offer a satisfying trade-off between accuracy and required computation time, which significantly increases with an increasing number of inclusions in the cell. The 50-inclusion samples in the present paper offer improved precision, while remaining tractable using modern hardware.
To determine average responses of the composites, the following analyses were performed. About 100 random microstructures were generated using the procedure described above for each investigated value of f i , which were: f i = 20% , 30% , 40% , and 50%. Out of each 100, two RVEs were selected, the first one having a microstructure with a very low value of λ ( Fig. 4 a) and the second one with a very high value of λ ( Fig. 4 b). Each sample was subjected to three isochoric tension tests, along each of the coordinate axes X, Y, and Z, under micro-periodic displacement boundary conditions. When extension was performed in the X direction, the global deformation tensor E of a sample was increased proportionally from 0 to E 11 = 0 . 030 , E 22 = E 33 = −0 . 015 , E 12 = E 13 = E 23 = 0 , and analogously for tests in the Y and Z directions. The FE analyses were performed in the AceFEM environment ( Korelc & Wriggers, 2016 ). To enforce the prescribed periodic boundary conditions resulting in the macroscopic (overall) strain E of the RVE, the displacement of every node of the FEM mesh lying on the surface of the cube and that of its periodic image on the opposite face of the cube were bound by the relation: where u A and u B are the displacements, while x A and x B are the initial positions, of the node ( A ) and its image ( B ), respectively. As each computation proceeded, volume-averaged stresses σ ij and volume-averaged deformations ε ij were recorded separately for each phase. On a desktop PC with the Intel i9-7900X processor and 64GB RAM, a single computation for an RVE with f i = 30% took between 6 hours for two-phase microstructures (see Fig. 3 ) and around 72 hours for three-phase microstructures (see Fig. 12 a). In the FEM analyses, microstructures were meshed with tetrahedral, 10 node, second-order finite elements, as in Pierard et al. (2007) , Czarnota et al. (2015) . The number of resulting nodes varied for different microstructures from 10 0 0 0 0 to 30 0 0 0 0, depending on the value of f i , the particular placement of inclusions, and especially on the presence of the interphase. The results of computations for a set of low-λ microstructures are shown in Fig. 4 a, and those for high-λ ones are shown in Fig. 4 b. In Fig. 4 , the Huber-von Mises stress is presented as a function of the strain tensor component in a given direction of elongation. For the inclusion volume fractions f i = 20% and 30%, the composite response is almost isotropic for both low and high λ values. The results for randomly generated composites presented in Section 4.2 concern MMCs with f i = 30% , for which the inclusions' parameters have a big impact on the effective microstructural response; still, 50 inclusions seem to be sufficient to have numerically isotropic behaviour of the composite.

Results
In this section, the basic two-pattern setup of the MRP model is validated for two-or three-phase isotropic particulate metal matrix composites against results of numerical homogenization performed on periodic unit cells and volume elements with random distributions of inclusions.

Regular arrangements of inclusions
In the numerical analysis, unit volume elements with regular arrangements of particles ( Fig. 5 a) were subjected to periodic boundary conditions (16) with a uniaxial strain tensor E : in which the unit vector v with components ( v X , v Y , v Z ) specifies the direction of stretching and ε is a scalar. For example, in the case of extension in the Z direction (along a cell edge) the components of v are (0,0,1), while for extension along a diagonal of the cube, they are 1 / √ 3 (1 , 1 , 1) . Representations of the overall strain tensor (17) in the basis coaxial with the unit cube's edges take the form respectively. In order to investigate the influence of the packing parameter itself, the two-phase metal matrix composite ( Table 1 ) will be used. The analyzed unit volume elements with regular arrangements of particles of the types Regular Cubic (RC), Body-Centered Cubic (BCC), and Face-Centered Cubic (FCC) are presented in Fig. 5 a. The notation RC, BCC and FCC is used because the periodic arrangements of spherical inclusions of equal radius R i in these cells follow the distributions of atoms in the RC, BCC and FCC crystal lattices Fig. 5 . The relations of the packing parameter λ / (2 R i ) to f i , valid for the respective regular arrangements, are as follows:  (17) ). The viewing direction is (1,1,1); due to the symmetry of the results, only single octants are presented. and the maximum volume fractions of non-overlaping inclusions, specified by the condition λ = 0 , are: The elastic-plastic response in terms of the σ Mises ( ε) curve obtained by FEM analysis of the assumed unit cells depends on the direction of tension, as seen in Fig. 5 b where the value of σ Mises at ε = 3% is shown. Unit cells exhibit extreme values of the equivalent Huber-von Mises stress σ Mises in tensile tests along: a direction corresponding to one of the cell edges (the Cartesian axes, (1,0,0), (0,1,0), and (0,0,1)), and a diagonal direction of the cube, e.g., (1,1,1) (see Fig. 5 b). Therefore, in Fig. 6 presenting full stress-strain curves, it was decided to show FE results for the respective unit cells in both (1,1,1) and (0,0,1) directions. It is worth noting that these are the directions of the extreme values of the directional Young's modulus determined analytically in Ostrowska-Maciejewska and Rychlewski (2001) , Cazzani and Rovati (2003) for materials with cubic symmetry. The Huber-von Mises equivalent stress σ Mises for ε = 3% , in the metal matrix composite undergoing a macroscopic tension test in a given direction v ( Eq. (17) ), is listed in Table 2 for the RC, BCC, and FCC unit cells with f i = 40% . Additionally, the mean value of σ Mises , calculated as the weighted average of results for all directions v over a unit sphere (see Fig. 5 b), is included in this table. It reaches its highest value for RC and the smallest for FCC. Note also that the difference between the minimal and maximal value is largest for the RC cell and smallest for the FCC cell. This difference can be treated as some measure of the anisotropy degree of the composite elastic-plastic response.
In Fig. 6 , estimations of the elastic-plastic response of metal matrix composites in uniaxial strain tests ( Eq. (17) ) are presented. The results of numerical homogenization are compared with estimations obtained from micromechanical models. In the analytical models under consideration, the overall response is independent of the extension direction and, in the case of the MRP approach, estimates are obtained as though for isotropic, random inclusion distributions for which the packing parameter follows the relations specified by Eq. (19) . The composite response is represented by the change in the equivalent Huber-von Mises stress with respect to the deformation ε in the direction of extension v . The results of micromechanical models shown in Fig. 6 were obtained using tangent linearization, because those for secant linearization differ significantly from the results of FEM. In Fig. 6 a-d, the influence of the inclusion volume fraction is shown, with the values of f i equal to 10%, 20%, 30%, and 40%. It is seen that with an increase in the volume fraction of reinforcements, the gap between SC and GSC estimates widens. The difference between the equivalent Huber-von Mises stress for f i = 10% obtained using the SC scheme or its generalization GSC is about 5 MPa at the strain ε = 3% , whereas for f i = 30% the difference between σ SC Mises and σ GSC Mises is around 65 MPa.
With an increase in the reinforcement volume fraction from 10% ( Fig. 6 a) to 40% ( Fig. 6 d), the difference in σ Mises obtained by numerical homogenization for directions (0,0,1) and (1,1,1) increases. Based on the analysis of the results presented in Fig. 5 b, it can be concluded that the stress-strain curves obtained for the other directions are between the two curves shown. In general, the results of numerical homogenization for all unit cells are below the results of micromechanical models with the exception of σ Mises ( ε) of the RC system for f i = 30% and 40% in the direction (0,0,1). For ceramic volume fractions f i = 10% and 20% , all FE results are below the micromechanical estimates ( Fig. 6 a and b). For f i = 30% , the FEM result for RC under extension along the coordinate axis (0,0,1) is similar to the generalized self-consistent estimation ( Fig. 6 c).
For f i = 40% , the estimation of σ Mises ( ε) obtained by the MRP scheme, for an isotropic distribution of inclusions and an identical packing parameter as in the corresponding unit cell, is consistent with the FEM response for the RC system under extension along one of the coordinate axes, i.e., in the direction of the maximum equivalent Huber-von Mises stress for this unit cell ( Fig. 5 a). However, the numerical results for the FCC and BCC systems are far below the micromechanical estimations. In Fig. 6 , in most cases, the elastic-plastic response of the metal matrix composite in the plastic flow range obtained with the GSC scheme is stiffer than the results of numerical homogenization. This is the opposite of what was found for elastic properties by Majewski et al. (2017) , for which, at any volume fraction of the inclusion phase, the elastic modules estimated with the GSC scheme are lower than those obtained using numerical homogenization.
The described relation between numerical and analytical estimates changes when an inclusion volume fraction above 40% is considered. In Fig. 7 , the influence of the inclusion volume fraction f i on the Huber-von Mises equivalent stress σ Mises is presented. As it was mentioned above, with an increase in f i the difference between the predictions of the GSC and SC schemes increases as it is shown in Fig. 7 a. In Fig. 7 b, the Huber-von Mises equivalent stress σ Mises for the strain ε = 3% is presented as a function of the inclusion volume fraction. Although for the volume fractions analyzed in Fig. 6 the analytical schemes overpredict the composite response, compared with the FE analysis, Fig. 7 b shows that for f i approaching the maximum values specified by Eq. (20) the FE calculations predict a stiffer response than the GSC scheme and the MRP approach, tending to the corresponding SC estimate. This is most readily seen when extension proceeds along the most densely packed direction in the unit cell. Size effect. Up till now, two-phase metal matrix composites have been studied. To investigate the size effect, an interphase (abbreviation 'int') is introduced as the third phase in the metal matrix composite, while the material parameters of the inclusion and matrix phases remain the same as in Table 1 . The interphase is a coating around the spherical inclusions, and its thickness t int is defined by the size parameter t int / R i . It is worth emphasizing that the inclusion volume fraction is kept constant, namely f i = 40% , so that when the interphase thickness increases, the volume fraction of the matrix decreases. The material parameters of the interphase are the same as those of the metal matrix with the exception of the yield strength, defined by the ratio Y int 0 /Y m 0 , which varies between 10 +12 and 10 −12 . For Y int 0 /Y m 0 < 1 , the geometry and material properties of the three-phase composite correspond to the problem of a degradation region between the reinforcements and the matrix. The reverse situation, Y int 0 /Y m 0 > 1 , corresponds to the case of a matrix reinforced in the close vicinity of inclusions. Note that for Y int 0 /Y m 0 = 1 the two-phase metal matrix composite case studied before is retrieved only when an additional layer is not introduced in the first pattern. Such a model (i.e. without an additional layer) is denoted in Fig. 8 as a 2-phase MMC.
The effects of the interphase quality and inclusion size in the metal matrix composite with the RC arrangement of inclusions are shown in Fig. 8 . The numerical results are obtained for a uniaxial strain test along the axis (0,0,1) ( Eq. (18) ). The results for a two-phase composite are plotted in Fig. 8 in green colour. For the assumed geometry of the two-phase composite, the MRP and FE estimates of the effective response are in good agreement. They are marked with a solid line (MRP) and discrete points (FEM), respectively. Additionally, in Fig. 8 a the local stress response in the interphase and two matrix phases is presented. The relative interphase yield strength Y int 0 /Y m 0 is assumed to be equal to 2 (blue curves) and 1/2 (yellow curves) in Fig. 8 a, and in Fig. 8 b to be between 10 +12 and 10 −12 , according to the legends. The size parameter is taken as constant, namely t int /R i = 1 / 20 in Fig. 8 a, while it varies from 10 −3 to 10 −1 (graph's X axis) in Fig. 8 b. Not surprisingly, as seen in Fig. 8 a, the presence of the interphase with enhanced properties Y int 0 /Y m 0 = 2 strengthens the macroscopic response of the composite and reversely, for the weakened interphase Y int 0 /Y m 0 = 1 / 2 , both the MRP and FEM results are below the estimates for the two-phase composite. It is worth noting that in the MRP model the influence of the interphase on the effective com posite response results not only from a change in the volume fractions of constituents, as in the SC scheme, but also from the interphase geometry that is modeled as a coating around the inclusions ( Fig. 12 a).
The influence of the interphase on σ Mises in the MMC estimated by the SC scheme is about half that predicted by the MRP approach, which is close to FEM outcomes ( Fig. 8 a). The micromechanical and numerical predictions of the macroscopic response of the three-phase composite are still in good agreement. However, differences can be noticed in the predictions of per-phase behaviour. It is worth mentioning that for all data in Fig. 8 a the stress is a function of the scalar ε of the macroscopic strain tensor E ( ε) ( Eq. (17) ). The metal matrix response obtained by numerical simulation is below the estimates obtained by the analytical MRP model, which distinguishes the part of the matrix around the composite spheres (the 4-GSC pattern) from the remaining matrix (the SC pattern). In Fig. 8 a, they are shown as the dash-dotted and dotted lines, respectively. The influence of the interphase yield strength on the matrix response in the 4-GSC pattern is clearly visible.   10. Contour maps of the Huber-von Mises equivalent stress for the isochoric tension test at the strain ε = 3% for the MMC composite with 50 spherical inclusions and f i = 30% ( Fig. 3 ). The packing parameter λ / (2 R i ) is equal to: (a) 0.039, (b) 0.164, (c) 0.260. Only the FE mesh of the matrix phase is shown.
A similar relation between numerical and analytical results is obtained for the interphase response itself: again, the MRP estimate is above the FEM results.
In Fig. 8 b, the size effect is explicitly presented. The horizontal axis t int / R i is on a logarithmic scale. Different ratios Y int 0 /Y m 0 are considered, as is indicated in the plot legend. Both approaches: analytical and numerical, show similar impact of the size parameter. When the ratio t int / R i tends to small values, the effective behaviour of the three-phase composite approaches the one of the two-phase material. The rate of decline of the MRP estimates is higher than that for the FEM results. Additionally, for t int / R i near the highest considered value 0.1, the differences between the FEM and MRP outcomes are significant.
In Fig. 9 , the influence of the yield strength ratio Y int 0 /Y m 0 is presented. Both the numerical or analytical results are obtained for the RC unit cell with different size parameters t int / R i . Due to numerical limitations of finite element analysis, only three size parameters t int / R i : 0.080, 0.050 and 0.015 were considered in the simulations. With an increase in the size parameter t int / R i , the difference between the results of numerical homogenization and micromechanical modelling also increases. Nevertheless, the results obtained with both methods show a similar character of the dependence of σ Mises on Y int 0 /Y m 0 . When the yield stress ratio tends to the limit values (in Fig. 9 , Y int 0 /Y m 0 on the X axis ranges from 10 −3 to 10 3 ), the σ Mises curves plateau and there is no further influence of Y int 0 /Y m 0 on the effective response of the MMC. In Fig. 12 of Majewski et al. (2017) , the influence of the interphase Young's modulus on the effective shear modulus of the composite is presented, in a manner analogous to the results in Fig. 9 of the present work. In the elastic range ( Majewski et al., 2017 ), if the interphase properties tend to 0 or infinity, the inclusion parameters lose their effect on the properties of the composite -the response of a three-phase composite tends to the response of a two-phase composite with voids or rigid inclusions. It is similar in the elastic-plastic range when the yield stress ratio exceeds a certain limit for a given interphase thickness; e.g., for t int /R i = 0 . 090 the limit ratio is about 10 -further changes in Y int 0 /Y m 0 do not affect σ Mises of the MMC because the interphase stays elastic in the process.

Random arrangements of inclusions
The procedure for generating RVEs with random arrangements of inclusions is described in Section 3 . The following results of numerical simulations were obtained for Metal Matrix Composites (MMC) reinforced with 50 randomly placed spherical inclusions with a volume fraction f i = 30% (see Fig. 3 ). The corresponding material parameters are provided in Table 1 .
Packing effect . To study the packing effect of inclusions on the effective elastic-plastic response of the composites, 100 different RVEs were generated. With their use, numerical simulations were performed of the isochoric tension test. Isochoric tension tests along the coordinate axes X, Y or Z are specified by the following strain tensors E , in the basis coaxial with the RVE edges: In the FE analyses, only 81 of the generated RVEs managed to reach the strain ε = 3% in the direction of stretching without numerical errors.
In Fig. 10 , contour maps of the Huber-von Mises equivalent stress for ε = 3% in Eq. (21)  These three microstructures are also presented in Fig. 3 .
For the smallest value of the packing parameter, namely λ / (2 R i ) = 0 . 039 ( Fig. 10 a), the maximum value of the Hubervon Mises equivalent stress in the MMC is higher -red colour according to the legend -than for the larger values of λ / (2 R i ) : 0.164 ( Fig. 10 b) and 0.260 ( Fig. 10 c). Let us recall that the inclusions are assumed to be elastic, whereas the matrix phase is considered as elastic-plastic ( Table 1 ). Due to the infinite contrast between yield stresses of the composite phases, extreme values of σ Mises are observed within these particles that cluster, i.e. when the distance to the nearest neighbour is small. This is clearly visible in Fig. 10 a in the upper part of the RVE, where three inclusions with the highest σ Mises values are located (dark red colours). By contrast, in the mean-field micromechanical models homogeneous mechanical fields are assumed in the dispersed phase. The examples of microstructures in Fig. 10 demonstrate that the stress fields in the inclusions are closer to homogeneous for larger values of the packing parameter. Fig. 11 a presents stress-strain curves for these three microstructures. The style of lines varies with the packing parameter λ / (2 R i ) . In Fig. 11 , estimations of the effective equivalent Huber-von Mises stress are shown, obtained by different micromechanical models: the SC, GSC and MT schemes, and the MRP approach using the packing parameters of the microstructures in Fig. 3 . In the initial stages of deformation ( ε < 0.5%), the effective equivalent Huber-von Mises stress determined by FEM is below the GSC estimates and above the MT estimates. In the further stages of deformation, the results of numerical homogenization are between the estimates of the self-consistent scheme and its generalization. The values of σ Mises for the analyzed structures are listed in the table inset in Fig. 11 a for both MRP and FEM approaches. Despite differences in the σ Mises results between the analytical and numerical approaches for the initial stages of deformation ( ε < 0.5%), the results for the final deformation level ( ε = 3 . 0% ) are in good agreement.
In Fig. 11 b, the effect of inclusion packing measured by λ / (2 R i ) on the overall composite response σ Mises at the strain ε = 3% in isochoric tension ( Eq. (21) ) is demonstrated. Analytical results for micromechanical models: SC, GSC, MT and MRP were obtained using tangent linearization. The outcomes of numerical homogenization are marked with points. In order to reveal the general character of the impact of the packing parameter λ / (2 R i ) on the effective material response, Fig. 11 b presents the results for 81 RVEs. For most of them, the isochoric tension test in two randomly selected directions: X, Y or Z ( Eq. (21) ) was simulated. The FE results for the RVEs allow one to estimate the impact of the packing parameter on the elastic-plastic response of the metal matrix composites. Because the SC and GSC models do not depend on the packing parameter λ / (2 R i ) , they result in horizontal lines in Fig. 11 b. When the packing parameter tends to 0, the MRP model tends to the SC scheme, whereas when the packing parameter tends to its maximum value (for the inclusion volume fraction of 30% it is about 0.49), the MRP approach tends to the GSC scheme. In spite of the significant scatter of the FEM results, they demonstrate a similar relationship: with a decrease in the packing parameter, on average, an increase in the overall equivalent stress is observed, and vice versa, for a larger value of λ / (2 R i ) , a lower effective stress σ Mises is observed. For larger values of the packing parameter the spread of the FEM results is smaller than for smaller values of the packing parameter. Although the MRP curve in Fig. 11 b differs from the FEM results in terms of σ Mises values, it can be concluded that the effect of inclusion packing on the effective composite response as predicted by the MRP model is qualitatively similar to that resulting from the performed FE analyses.
Size effect. The inclusion size parameter is determined by the ratio t int / R i . In Majewski et al. (2017) , the impact of the inclusion size on the effective shear modulus was investigated. Unfortunately, carrying out similar tests for the elastic-plastic range is much more difficult. The addition of an interphase, even of a relatively large thickness t int /R i = 0 . 20 , increases the   number of nodes several times, extending the time of a single numerical simulation from several hours, when without an interphase, to about three days. Therefore, it has been decided to present the results of numerical simulations only for a selected RVE.
The respective microstructure of the two-phase composite is shown in Fig. 3 c (the packing parameter equals 0.260). In Fig. 12 a, contour maps of the Huber-von Mises equivalent stress are presented for the isochoric tension test at the strain ε = 3% after adding the interphase coating. Only the FE mesh of the interphase is shown. The legend for the contour map is the same as in Fig. 10 . In the first case ( Fig. 12 a top), the yield stress of the interphase is twice as high as that of the matrix: Y int 0 /Y m 0 = 2 , while in the second case ( Fig. 12 a bottom) an inverse relation is prescribed: Y int 0 /Y m 0 = 1 / 2 . The presence of a hardened zone around inclusions increases the stress in the inclusions compared with the case of the two-phase composite ( Fig. 10 c). In turn, the introduction of a weakened zone between the inclusion phase and the matrix reduces the stress in the inclusions. The trend follows the one observed for a unit cell with a regular arrangement of particles studied in the previous subsection.
In Fig. 12 b, the influence of the interphase (the size parameter is 0.20) on the effective response of the MMC is presented. The green color indicates the macroscopic response of the two-phase composite, the blue color denotes the MMC with a reinforcement in the form of interphase coatings around inclusions ( Y int 0 /Y m 0 = 2 ), while the yellow color corresponds to the case of a weakened interphase ( Y int 0 /Y m 0 = 1 / 2 ). As before (see Fig. 8 ), the presence of the additional phase affects the σ Mises ( ε) curve. The results of MRP modelling and numerical simulations are in good agreement for the final deformation ε = 3% , while at the initial deformation stage the difference between micromechanical estimates and numerical homogenization is significant, especially in the case of the hardened coating. Fig. 12 c presents the σ Mises ( ε) curves for individual phases of the three-phase composite. Good agreement between the MRP and FEM predictions for the interphase is found. However, the difference between MRP and FEM in the estimated state of stress in the matrix is significant, similarly to the RC unit cell ( Fig. 8 a), which probably results from the uniform field assumption in the MRP approach. Due to the observed analogy between the results for regular and random arrangements of inclusions in the selected example it is expected that variation of the size parameter leads to a change in the yield stress analogous to the one shown in Fig. 9 .

Conclusions
The influence of the particle size and packing on the overall elastic-plastic properties of metal-matrix composites has been analyzed. The mean-field Morphologically Representative Pattern (MRP) approach, originally developed to address linear material behaviour, has been adapted to account for non-linearity. Within the proposed MRP formulation, a particulate composite microstructure has been represented by two patterns: one with a spherical inclusion surrounded by a layer having a thickness defined by the mean minimum distance between particles in the representative volume, and the second one carrying the remaining pure matrix material. Both patterns are embedded in a homogenized effective medium. The mean minimum distance parameter describes the packing (clustering) of the particles. To take into account the size effects, an interphase surrounding the inclusion in the first pattern has been introduced. Its thickness plays the role of the size parameter.
The incremental linearization, proposed by Hill, is applied in order to take advantage of the mean-field models for elasticplastic materials. Two variants of the linearization procedure: tangent and secant have been considered, depending on the definition of the current stiffness tensor. The metal matrix is assumed to be a ductile material with linear elasticity and the Huber-von Mises yield function, with the associated flow rule, whereas the inclusion phase is assumed to be a linearly elastic material with parameters relevant to a ceramic. As it is commonly observed for mean-field frameworks extended to elasto-plasticity ( Kursa et al., 2018 ), the MRP approach combined with secant linearization predicts a stiffer response of an MMC than when it is combined with tangent linearization.
The validity of the analytical mean-field estimates, and in particular of the MRP approach, with respect to computational homogenization results has been analyzed. The numerical simulations have been performed using the finite element (FE) method. Two types of composite volume elements have been subjected to numerical homogenization: periodic unit cells of cubic crystal-type arrangements and representative volume elements with a random distribution of particles. In the first case, the elastic-plastic MMC response depends on the loading direction. In the uniaxial strain process in different directions, the highest average value of the Huber-von Mises equivalent stress is obtained for the regular cubic placement of particles and the smallest one for the FCC cell. At the same time, the largest and the smallest anisotropy degree (assessed in terms of the difference between the minimum and maximum values obtained for a given cell) have been respectively observed for these two unit cells.
It is noticed that with an increase in the volume fraction of ceramic particles the gap between the SC and GSC schemes increases. In general, the results of numerical homogenization for a given unit cell are below the results of micromechanical models with the exception of the volume fraction of reinforcements above 40%.
In order to reveal the general character of the impact of the packing parameter on the effective material response, over 80 RVEs with 50 randomly placed particles have been generated and the isochoric tension test has been performed. Locally, due to the infinite contrast between yield stresses of the composite phases, extreme values of σ Mises for a given strain level have been observed in FE calculations within these particles that cluster, i.e. when the distance to the nearest neighbour is small. As concerns the overall response, in spite of the significant spread of the FEM results they show a similar effect of packing to the one predicted by the MRP model: with a decrease in the packing parameter, on average, an increase in the overall equivalent stress is observed, and inversely, for a larger value of the packing parameter, a lower effective stress is observed. For larger packing parameters, the scatter of the FEM results is smaller than for small values of the packing parameter. Although the MRP predictions differ from the FEM results in the value of the Huber-von Mises equivalent stress, it can be concluded that the effect of inclusion packing on the effective composite response according to the MRP model is qualitatively the same as in the performed FE analyses. It should be noted that the standard Mori-Tanaka, SC and GSC schemes cannot account for the effect of reinforcement packing in particulate composites.
To investigate the size effect, an interphase has been introduced as a third phase in the metal-matrix composite in the form of a layer around the spherical inclusions. The inclusion volume fraction has been kept constant. The material parameters of the interphase are the same as of the metal matrix, with the exception of the initial yield stress. Not surprisingly, the presence of the interphase with enhanced properties ( Y int 0 /Y m 0 > 1 ) strengthens the macroscopic response of the composite and reversely, for the weakened interphase ( Y int 0 /Y m 0 < 1 ), both the MRP and FEM results are below the estimates for the two-phase composite. The influence of the interphase on the effective composite response is not only due to the change in the volume fractions of constituents but also due to the interphase geometry. The numerical and micromechanical modelling of the overall response of the three-phase composite are in good agreement; however, differences are noticed in the predictions of per-phase behaviour.
The performed study indicates that the MRP approach can be used as an effective predictive alternative to computational homogenization, not only in the case of linear elasticity, as it was shown earlier ( Majewski et al., 2017 ), but also in the case of non-linear elastic-plastic composite behaviour.