Yield Surfaces of Material Composed of Porous and Heterogeneous Microstructures Considering Phase Debonding

This work deals with numerical simulation of the mechanical behavior of materials composed of heterogeneous ductile microstructures using a multi-scale approach considering plasticity processes and phase debonding. Due to few studies about yield surfaces of metal matrix composites (MMC) with weak interfaces presented in the literature, the major goal of this work is to propose yield surfaces for metal matrix composites reinforced by rigid inclusions. The yield surfaces are obtained for Representative Volume Elements (RVEs) of materials presenting perfectly bonded inclusions and phase debonding in the interface zone. The matrix is considered an ideally plastic material governed by von Mises model, whereas the interface zone is modeled by means contact and fracture constitutive models incorporated in a proposed finite element. Also, RVEs containing different distributions and volume ratios of voids are analyzed. Considering the phase debonding, for compressive loadings the RVE behaves like RVE with perfectly bonded inclusions whereas for tension loadings the RVE presents a behavior quite similar to the one with voids. On the other hand, the concentration of voids in the RVE decreases its mechanical strength.


INTRODUCTION
The knowledge of the mechanical properties of the materials is an important aspect for the Materials Science. In this context, computational modeling is an attractive alternative, since it is possible to obtain results using numerical simulations replacing several costly laboratory tests. Other aspect is related to the capabilities increasing of the computational devices allowing the development of the constitutive models and numerical techniques. Due to the similar characteristics of the mechanical behavior at macroscopic level, many constitutive models based on continuum mechanics and thermodynamics of solids have been proposed for limited applications in several cases, Peric et al. (2010), Nguyen et al. (2011) and Pituba et al. (2016). This limitation is evident when dealing with materials composed of heterogeneous microstructures resulting in constitutive models with complex formulation and many parameters to identify (Karamnedjad et al., 2013). However, deformation and rupture processes take place at microscopic level. In this context and taking into account the advances on the computational mechanics, many numerical techniques and constitutive models have been proposed to describe the mechanical behavior of heterogeneous materials at microscopic level considering fracture and plasticity processes, for example, Pituba and Souza Neto (2015), Azizi (2012), Fernandes et al. (2015a), Fernandes et al. (2015b), Santos et al. (2016), Nguyen et al. (2011), Cavalcante et al. (2011, Needleman and Tvergaard (1987), Tvergaard and Needleman (2006)among others.
In this work, a finite-element procedure proposed by Peric et al. (2010) and extended by Pituba et al. (2016) to take into account the fracture process in the EVR within a purely kinematical multiscale framework is used to simulate the mechanical behavior of metal matrix composites considering the influence of the phase debonding. For Giusti et al. (2009), the need for more precise responses has improved the studies on the lower scales leading to the development of multiscale theories, which use informations from different scales to obtain the material response. According to Miehe (2003) and Reis and Pires (2013), the use of homogenization techniques for modeling heterogeneous media is an interesting alternative and important subject of the research field. The initial focus is to understand the physical mechanisms that are developed in the microstructure of the material and, later, this information is used to obtain results about complex macroscopic behaviors. At each point of the macroscopic level, a local analysis is carried out through the modeling of a Representative Volume Element (RVE), which is responsible for describing the material microstructure. The micro-to-macro transition is based on the Hill-Mandel principle, which establishes the existence of the energy equivalence between a point on the macro-continuum and the RVE. In this context, studies on quasi-brittle, composites and porous ductile materials are included, mainly in which, heterogeneities significantly influence the behavior of the material. Finally, as the main advantage of multi-scale theories, it is possible to highlight the more detailed definition of the material microstructure and, consequently, a more detailed analysis of its behavior.
On the other hand, the composite materials are formed by the combination of two or more constituents with different properties resulting in a material with specific characteristics. Thereby, it is expected to obtain a material with improved properties when compared to the individual constituents. Accordingly to Lewis (2014), the Boeing 787 is composed of more than 50% of composite materials giving us an example of the importance of using this class of materials. The mechanical behavior of composite materials is affected by several parameters, including the size and proportion of the inclusions, as well as their distribution in the matrix. In this work, special attention is given to another very relevant aspect related to the matrix-inclusion interface: the phase debonding. Accordingly Tian et al. (2010), the fracture process in the composite materials involves the nucleation of voids at the interface between the constituents and, later, the growth and coalescence of these voids in the matrix zone. A weak interface results in low stiffness and strength. This decreasing process of the mechanical properties occurs because the weak interfaces contribute significantly to the fracture process of the material through the matrix-inclusion decohesion. It is noteworthy that several studies have been developed to model the decohesion process at the interface zone in composite materials, such as Ghosh et al. (2000), Ghassemieh (2002), Chen et al. (2003), Sun et al. (2003), Papakaliatakis and Karalekas (2004), Segurado and LLorca (2005), Zhang and Chen (2012), Abadi (2012) and Azizi (2012).
However, a research field that still requires studies includes the proposal and study of beginning and evolution criteria of the plasticity process in heterogeneous media presenting ductile behavior. This assertion is due to many constitutive models have been proposed in the past for materials considering some requirements, such as: homogeneity, continuum medium, isotropy and incompressibility. For some cases, these requirements do not lead to realistic results when dealing with ductile materials with inclusions and voids. For this reason, some works have been recently developed in order to describe yield surfaces for heterogeneous materials. For ductile media with small voids proportions, Gurson (1977) has proposed a yield surface within an analytical modeling. This pioneering work has been the reference for others works dealing with yield surfaces for heterogeneous media. For example, Tvergaard (1981) has proposed an improvement on Gurson's model to deal with shear band instabilities. Also, Needleman and co-authors have studied ductile fracture processes in microstructures with random distributions of void nucleation in 2D and 3D analyses (Needleman and Tvergaard (1987) and Tvergaard and Needleman (2006)). Brünig, Voyiadjis and others have used homogenization techniques for modeling damage processes in composite materials (Voyiadjis and Kattan (1993), Kattan and Voyiadjis (2003), Brünig et al. (2011) andBrünig et al. (2014)). Another work dealing with porous metals has been developed by Giusti et al. (2009), which used a multiscale approach to obtain results considering different void ratios, being compared with ones presented by Gurson (1977).
Regarding to ductile materials with inclusions, the work of Gărăjeu and Suquet (1997) deserves to be highlighted. The authors have presented results for viscoelastic and plastic matrix containing ideally rigid particles using a semi-analytical modeling considering a perfect bonding between matrix and inclusions. This situation does not properly represent the mechanical behavior of composite materials presenting a weak interface. Recently, the subject has been taken up by the work of Somer et al. (2015), which estimates yield surfaces for metals with inclusions using a multiscale approach considering phase debonding. The authors have analyzed RVEs containing an inclusion and the interface has been modeled by a Coulomb type fiction law. In this context, this work intends to contribute to the knowledge about yield surface for metallic materials with inclusions perfectly bonded and the influence of the weak interface in the yield surface using simple constitutive models at microscopic level within a purely kinematical multiscale framework. For this propose, the weak interface is modeled by a recently proposed contact and fracture model in Pituba et al. (2016). This model represents the decohesion process in matrix/inclusion interface leading to a proper representation about the homogenized macroscopic behavior of the material. Besides, this work deals with yield surface for materials containing different distributions and volume fractions of voids.
The paper is divided in five sections, as follows: a brief description on multiscale framework is presented in Section 2 as well as the incorporation of the fracture process in the formulation. In Section 3 is discussed the constitutive models employed to represent the dissipative phenomena at Latin American Journal of Solids and Structures 14 (2017) 1387-1415 microscopic level. Section 4 presents numerical analyses and yield surfaces proposed for metal matrix composites considering inclusions perfectly bonded and the influence of the weak interface consideration. Moreover, yield surfaces for porous ductile materials represented by RVEs containing voids are proposed. Finally, some discussions and possible extensions are presented in Section 5.

MULTISCALE FORMULATION CONSIDERING FRACTURE PROCESS. OVERVIEW
In the context of the multiscale approach, each point x of the macroscale W is considered as a RVE representing the microscale of the material m W . Each point over the volume of the RVE is characterized by coordinate point y and the boundary of the RVE is denoted by m ¶W . Besides, the length of the microscale has to be smaller than the length of the macroscale ( ) l l m  . However, it is important to characterize the proportionality between the material phases over the RVE which composes the microstructure of the material. The macroscopic properties of each point of the macroscale are obtained by means homogenization techniques applied on the RVE volume. An example of multiscale analyses is shown in Figure 1.
The RVE is considered as a continuum medium, so that the stress concept is valid at microscale. The macroscopic quantities for strain ( ) Equations (1) and (2) represent the macroscopic or homogenized values for strain and stress, as a microscopic filed have been transformed into a macroscopic quantity by means of a homogenization technique.
On the other hand, by the homogenization process it can also obtain the homogenized constitutive tangent modulus ep C , as follows:  (2) and (3). Besides, the microscopic stress can be written in terms of the microscopic strain, as follows: where y f is the constitutive functional. In this work, for the triangular elements defined in the matrix (see section 3), the constitutive functional y f is defined by the von Mises elasto-plastic criterion while for the rectangular elements defined on the interface between matrix and inclusions, the stress are computed taking into account the fracture and contact phenomena. Moreover, the microscopic strain m e can be written in terms of the microscopic displacement filed m u of the RVE, as follows: where S  is the symmetric gradient operator of the displacement field u . Without loss of generality, the microscopic displacement filed m u can be defined as the sum of three parts: being the first one a constant representing a rigid body motion coincident to the macroscopic displacement ( , ) t u x related to the point x , the second one is obtained from the macroscopic strain e as follows: which varies linearly with the coordinate y , and the last part is a displacement fluctuation field as the macroscopic displacement ( , ) t u x is a rigid body motion, it has no influence in the stress field at points y of the RVE, therefore it is not taken into account to obtain the solution of the equilibrium problem. Thus, in equation (8) ( , ) t u x has been disregarded, see Giusti (2009).
In Equation (8) the part ey varies linearly with y resulting from the multiplication of the macroscopic strain e of the RVE, which is constant, by the coordinates of the point y . In the case of having uniform microscopic displacement m e , the displacement fluctuation  m u is null. In the RVE the following relations for the microscopic strain m e and the microscopic strain fluctuation  m e have to be satisfied: Considering Equations (8) to (10) the microscopic strain can also be written as: After some manipulations (Fernandes et al., 2015a), Equation (11) can be written in terms of velocity, where a microscopic strain velocity is cinematically admissible if: where m n is the space of cinematically admissible displacements of the RVE. More details can be found in Fernandes et al. (2015a). As already mentioned, the microscale is represented by the RVE, being the FEM (Finite Element Method), the numerical method used to solve the RVE equilibrium problem. Therefore, the computation of displacements, internal forces, stress and constitutive tensors, for all finite elements, are obtained when the convergence of the equilibrium problem is achieved according to the adopted tolerance factor. But, in order to solve the RVE equilibrium problem, boundary conditions in terms of displacement fluctuations must be imposed to the RVE. Therefore, considering Equations (3), (9) and (11), the following Equation can be obtained to represent the equilibrium problem of the solid part of the RVE in terms of displacement fluctuation: Finally the formulation is completed by the appropriated choice of the space m n , i. e., with the choice of the kinematical restrictions to be imposed to the RVE. Thus, the microscopic equilibrium problem consists of, given the macroscopic strain tensor e , finding the field for iteration i+1, such that: where F is the force vector and K the tangent stiffness matrix of the RVE. After computing the (14), the next step is to obtain the displacement fluctuation field to be considered at iteration i+1 given by: Accordingly Pituba et al. (2016) On the other hand, the homogenized stress is computed from Equation (2), considering that the RVE is composed by voids and a solid part (matrix and rigid inclusions) The RVE equilibrium problem is completed with the choice of the kinematical restrictions to be imposed to the RVE, leading to different classes of multiscale models and consequently to different numerical results (Peric et al., 2010). In this work only periodic displacement fluctuations is considered. For that, each RVE side i + G whose normal direction is i n + , must correspond to an equal side

CONSTITUTIVE MODELING
In order to simulate the plastic strains presented in the microstructures of the ductile porous materials, the von Mises model with strain hardening has been used. Besides, for analysis of MMC, cohesive fracture and contact models have been applied to describe the phase debonding in the interface zone as well as a rigid and elastic behavior is adopted for the inclusion. The Finite Element Method (FEM) is used to model the RVE problem. Triangular finite elements have been employed in the discretization problem of the inclusions and matrix zones. Contact fracture finite elements recently developed in Pituba et al. (2016) have been used in the interface zone. Figure 2 shows the discretization scheme. The nucleation and evolution of the fracture process have been added to the previous microstructure formulation (item 2), leading to different homogenized stress and constitutive tensors. For the cohesive fracture model proposed in Pituba et al. (2016) and used in this work, the traction and opening displacement relation considers a gradual separation process of the fracture surfaces in order to avoid numerical instabilities due to the strong discontinuity. Besides, each Gauss point of the contact fracture finite element (Figure 3) contributes to the internal force evaluation by means of the traction vector computed either by the cohesive law (if a crack is opened at that Gauss point) or by the contact law (if a crack is closed at that Gauss point). In the coordinate local system of the contact fracture finite element (Figure 3), the s and n axes indicate the sliding and normal direction, respectively.

()
(2 ) ( 4 ) It is important to note that the contact fracture finite element is composed of two surfaces which are coincident in the undeformed configuration of the RVE. The contact fracture finite element is defined as an element with four nodes and its geometry is compatible with the two triangular finite elements used to model the matrix and inclusion zones. The formulation of the contact fracture finite element is presented in Pituba and Souza Neto (2015) and Pituba et al. (2016). Pituba et al. (2016) have proposed a cohesive fracture law in order to deal with damage process leading to the complete failure of microstructures in ductile media. The constitutive model describes the finite-deformation irreversible cohesive law. The cohesive free energy is given by: where, n d is the normal opening displacement due to mode I; S d is the sliding opening displacement due to mode II and q is the internal variable that describes the inelastic processes related to decohesion. It is possible to assume that the deformation due to sliding opening process is a scalar value independent of the sliding direction on the cohesive surface, thus S d = S d , therefore the behavior has an isotropic characteristic and the cohesive law is written introducing an effective opening displacement expressed by: The parameter b assumes different values (from 0 to 1) establishing a weight ratio between the sliding and normal directions. The f free energy potential depends of d , and the cohesive law is where, n is the unit normal to the cohesive surface; S d is the sliding opening vector located on the cohesive surface, t is the cohesive traction on the crack; t is a scalar effective traction.
On the other hand, the released cohesive energy in the microstructure of the material proposed in this work (Equation (20)) is given by: where the law for the scalar effective traction for the loading cases is obtained from Equation (23) as: For the scalar effective traction for the unloading cases is proposed a law considering an elastic behavior, i. e., without residual effective opening displacement as follows: Where e is the e-number, c s is the maximum tension cohesive normal traction and c d is a characteristic opening displacement that indicates a critical opening. Thus, b , c s and c d are parameters of the cohesive model. Besides, d  is the opening displacement rate. Also, there is a relation between the cohesive law and the critical energy released rate (GC) for crack propagation in the microstructure. Assuming the direction 1 as the direction on the fracture surface and towards to the its propagation, GC can be written as: where R is the cohesive zone length. The Equation (26) can also be defined as: For the cohesive law presented in this work, using Equation (24), the critical energy released rate is given by: Before the nucleation process, it is assumed the existence of stiffness between the surfaces of the future fracture situated between triangular finite elements. This stiffness is simulated by another parameter of the proposed model called penalty factor ( p l ). In a practical view, high values for this parameter are adopted in order to obtain an accurate approach. This procedure ensures that the future fracture be kept closed until the separation criterion is reached and, at the same time, guarantees the physical admissibility of the entire process. The penalty factor is, therefore, a stiffness imposed to the crack keeping it closed.
In general way, that strategy intends to create stiffness between the nodes of the embedded contact fracture finite elements in the interface zone in order not to allow penetration between the surfaces of the fracture. However, in tension regimes, this penalty factor effectively replaces the initial rigid part of the cohesive law for a linear response given by Equation (29). In order to detect the cohesive contact phenomenon, the concept of the opening displacement gap between the Gauss points of the contact fracture finite element is adopted.

RESULTS
In this section, the computational homogenization-based approach described in section 2 is used for numerical analyses of MMCs (metal matrix composites) and porous ductile microstructures in order to obtain yield surfaces and understanding the mechanical behavior of these kinds of materials. For the MMCs, the analyses are performed considering perfect bonding in the interface zone as well as the phase debonding to discuss the influence of this fracture process on the mechanical strength of MMCs  q obtained in each analysis have been used for the construction of yield surfaces for the materials.

Numerical Analysis of Metal Matrix Composites
This works intends to discuss the progressive collapse of heterogeneous materials composed of ductile microstructures reinforced by rigid inclusions and weaken by voids. Somer et al. (2015) has discussed this theme using a Coulomb type fiction law to model the interface. However, in this section, cohesive and contact laws are used with contact fracture finite element. In Pituba et al. (2016) is presented some advantages of the proposed modeling, such as: the proposal of a smoother stiffness reduction for ductile media, and the stiffness of the contact fracture finite element is computed taking into account individual contributions of each Gauss point, which can be either in an opened or closed condition, as described in section 3.
In order to observe the mechanical behavior of metal matrix composites, different RVEs submitted to different loading states have been analyzed. For the study, four different RVEs with their dimensions defined by LxL and thickness given by L/10 have been used: one RVE containing a concentrated void at the center of the RVE with 10% porosity; another RVE containing an inclusion, at the center of the RVE, perfectly bonded to the matrix with 10% of volume fraction; the same configuration of the last RVE, but considering the phase debonding with b = 0 (considering only normal opening of the crack surface) and b = 0.707 (considering normal and sliding phenomena over the crack surface).
For the RVE with void 1800 triangular finite elements have been used. For the RVE with an inclusion perfectly bonded, 1976 triangular finite elements have been used. Finally, for the two situations of the RVE containing an inclusion considering phase debonding, 1976 triangular finite elements have been used and 32 contact fracture finite elements have been employed to model the interface zone. In all numerical analyses, plane strain conditions in small strain regime have been considered using a tolerance factor of 10 -6 to check the convergence of the nonlinear procedure.
On the other hand, von Mises Model with perfect elasto-plasticity has been adopted to model the mechanical behavior of the matrix zone considering a yield stress s = To perform the analyses presented in this section, a macroscopic strain tensor have to be imposed to the RVE. In order to cover a range of possible loading states, the Equation (32) is adopted to describe the imposed macroscopic strain tensor: where a represents a load factor. Besides, in all numerical analyses the macroscopic strain tensor has been divided in 30 increments to perform the non-linear analysis (∑ N° Increments = 30, see Figures 4,5,6,8,9 and 10).
Therefore, the normal strains ( xx e and yy e ) and/or distortional strain ( xy g ) vary accordingly with adopted value for a . For the analyses focusing on behavior of the normalized effective stress ' q and pressure ' p , positive ( a = 1.0 , a = 0.8 and a = 0.6) and negative ( a = -1.0 , a = -0.8 and a = -0.6) values for a are used. Besides, the analysis with a = 0.0 has been performed, but for these cases the numerical convergence has not achieved due to the high values for the distortional macroscopic strain imposed to the RVE leading to early plastic collapse of the matrix zone.       For  = -0.6 Figure 11 presents the x-direction stress distribution over the RVE for negative values of load factor. The same observation reported for Figure 7 about the visualization of the stress distribution is valid here for a = -1.0 .
Also, more analyses were performed using others values for α considering positive and negative values in order to provide a range of results which are necessary to propose the yield surfaces related to obtained maximum values for ' q and ' p . The yield surfaces considering microstructures with void, bonded inclusion and debonded inclusion are presented in Figure 12. In order to compare the results obtained by the proposed modeling for porous ductile materials with Gurson's model, Figure  12 also shows the results expressed by Equation (33) proposed by Gurson (1977).
where f represents the void volume fraction (in this case, ) and ' p and ' q are given by Equations (30) and (31), respectively. The numerical results show that for compressive loadings the RVE considering phase debonding presents a similar mechanical behavior to the RVE with bonded inclusion. This conclusion is in agreement with Somer et al. (2015). This result was expected because in case of compressive loadings the contact phenomena play an important role in the mechanical behavior of the RVE keeping the possible fractures closed. Therefore, even considering phase debonding, the inclusion is capable to conferee an high rigidity to the RVE. Therefore, for the cases of RVE with bonded inclusion and RVE considering the phase debonding, a significative parcel of the stress state is supported by the reinforcement (rigid inclusion).
For expansive loading cases the RVE considering phase debonding presents a similar mechanical behavior to the RVE with centered void. This qualitative result is also in agreement with Somer et al. (2015). This similarity of mechanical behaviors is possible due to the gradual separation of the fracture surfaces leading to a complete debonding and failure of the microstructure as a RVE with void, i. e., the rigid inclusion does not play its role as reinforcement. It is important to note that for high values of distortional strains, difficulties to obtain the numerical responses of the RVEs considering phase debonding have been experienced, mainly for the case with b = 0 where the shear stresses have leaded to a massive yielding of the matrix zone. Besides, it is important to note that the proposed modeling presents results quite satisfactory in expansive loading cases for microstructures with phase debonding or voids when compared with Gurson's model. Moreover, for microstructures with voids in compressive loading cases, the proposed modeling also is agreement with Gurson's model. Finally, based on the results presented in this section some approximate equations have been proposed to deal with the different kinds of heterogeneous microstructures. Figure 13 shows the proposed equation for reinforced microstructures considering perfectly bonded inclusion. Figures 14  and 15 show the proposed equations for reinforced microstructures considering phase debonding ( b = 0 ) submitted to expansive and compressive loadings, respectively. Figures 16 and 17 are about the consideration with b = 0.707 . Figure 18 presents the proposed equation for porous ductile materials.

Numerical Analysis of Porous Ductile Materials
In this section, numerical analyses of porous ductile microstructures are presented considering the following models: i) RVE containing a concentrated void with 5% porosity where 1908 triangular finite elements have been used; ii) RVE containing randomly distributed voids with 5% porosity where 2117 triangular finite elements have been used; iii) RVE containing a concentrated void with 10% porosity where 1800 triangular finite elements are used; iv) RVE containing randomly distrib-q' = -0.00005932p' 6 -0.00233397p' 5 -0.03745418p' 4 -0.31277838p' 3 -1.43600419p' 2 -3. uted voids with 10% porosity where 1922 triangular finite elements have been used. The major goal here is to analyze the differences of the mechanical behavior of porous materials considering two factors: the value of the porosity and the distribution or concentration of voids.  In all analyses, Von Mises Model with perfect elasto-plasticity has been adopted to model the mechanical behavior of the matrix zone considering a yield stress y s = 240 MPa . Also, a Young´s where a represents a load factor. Now, RVEs containing porous ductile material are analyzed, therefore it was necessary to decrease the values for the imposed macroscopic strains in order to avoid the massive yielding of the matrix zone leading to the collapse of the RVE and, also, in order to obtain the convergence of the numerical procedure, mainly in the RVE with 10% porosity with distributed voids. Besides, in all numerical analyses the macroscopic strain tensor has been divided in 30 increments to perform the non-linear analysis (∑ N° Increments = 30,see Figures 19,20,21,22,24,25 and 26).
For the analyses focusing on behavior of the normalized effective stress ' q and pressure         Figure 27 presents the x-direction stress distribution over the RVE for negative values of load factor.
As explained in section 4.1, more analyses were performed using others values for a in order to propose the yield surfaces considering microstructures with randomly distributed voids and concentrated void presented in Figure 28.   In general way, it can be noted that the mechanical behavior of the RVEs with a concentrated void is similar to the RVE with randomly distributed voids related to component ' p for expansive and compressive load factors (see Figures 19,20,21,24,25 and 26). Although the results are quite similar, the RVEs with randomly distributed voids present higher values of ' p than the RVEs with concentrated void. As expected, the RVEs with 5% porosity have presented higher strength than ones with 10% porosity. For 0.0 a = , the numerical results were different and small values have been observed being considered nulls (see Figure 22) because in this case only macroscopic distortional strains are imposed leading to small values for ' p , mainly in RVEs with symmetric distributions of voids.
On the other hand, for component ' q , it can observed that the void distribution has influenced the numerical results. In general point of view, for the same void ratio, the values for ' q in RVEs with randomly distributed voids are bigger than for the cases of RVEs with concentrated voids. Therefore, the concentration of voids at center of the RVE has influenced the deviatoric part of the stress tensor leading to decreasing of values for component ' q . The RVEs with 5% porosity have presented values for ' q bigger than RVEs with 10% porosity in all cases, except in situations with a = 1.0 and a = -1.0 . In these last cases, the numerical results have presented small values for ' q due to only macroscopic normal strains have been imposed to the RVE.
Besides, Figure 28 shows that the increasing of the porosity implies in yield surfaces with smaller values. The increasing of the porosity leads to the loss of strength of the RVE and its early collapse. This conclusion is in agreement with Giusti et al. (2009). Moreover, the shape of the void distribution has influenced the obtained results. The concentration of voids has implied in yield surfaces with smaller values when compared to distributed voids situation. It is important to note that the yield surfaces have presented little different values for expansive and compressive loadings. It is possible to observe that expansive loadings present a more evident difference between surfaces when considering concentrated and distributed voids, in each situation (5% and 10% porosities).
For 5% porosity the following equation has been proposed based on the results presented in Figure 28: ).

CONCLUSIONS
In this work yield surfaces for porous ductile media and metal matrix composites considering phase debonding have been proposed. A computational homogenization-based approach considering frac-ture processes has been used to model the mechanical behavior of materials at microscopic level. The importance to consider the phase debonding in metal matrix materials has been addressed in Pituba et al. (2016). However, it is necessary to discuss the safety limits to deal with this kind of materials. Therefore, this work has contributed to understand the mechanical behavior of MMCs and porous ductile materials at microscopic level leading to properly proposed yield surfaces.
Specific results have been found in this work related to MMC. For instance, it is possible to note that for compressive loadings, RVE considering phase debonding has the same mechanical behavior of the RVE considering perfectly bonded inclusion. However, for expansive loadings, RVE considering phase debonding behaves closely to the RVE with void. Therefore, the consideration of phase debonding is important in collapse regimes when the material is excited by dominant tension regimes.
Another interesting discussion is about the porous ductile media. In general way, the numerical results have shown that the increasing of the void volume fraction leads to yield surfaces with small values. Besides, the distribution of voids in the RVE has a massive influence on the Von Mises effective stress. If the void volume fraction is kept constant, the RVEs with randomly distributed voids present high values for the Von Mises effective stress when compared to the RVEs with concentrated void. So, it is possible to observe that the concentration of voids in the RVE leads to the decreasing of its strength.
Finally, the computational homogenization-based approach has shown the capability to deal with complex macroscopic phenomena of mechanical behavior of heterogeneous materials using simple, but efficient constitutive models at microscopic level. This formulation is potentially applicable to fully coupled multiscale analysis of solids composed of heterogeneous materials.