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

Abstract

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.

Keywords:
RVE; cohesive zone model; homogenization; phase debonding; yield surface; metal matrix composites

1 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) and (Brünig et al. 2014)). Another work dealing with porous metals has been developed by (Giusti et al. 2009), which used a multi-scale 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 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.

2 MULTISCALE FORMULATION CONSIDERING FRACTURE PROCESS. OVERVIEW

In the context of the multiscale approach, each point x of the macroscale Ω is considered as a RVE representing the microscale of the material Ω_μ . Each point over the volume of the RVE is characterized by coordinate point y and the boundary of the RVE is denoted by ∂Ω_μ . Besides, the length of the microscale has to be smaller than the length of the macroscale (_{l ≫ lμ} ). 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.

Figure 1
Multiscale analysis scheme considering two scales.

The RVE is considered as a continuum medium, so that the stress concept is valid at microscale. The macroscopic quantities for strain ε (x, t ) and stress σ (x, t ) at a point x of the macro-continuum are defined as the volumetric average of their respective field _{^{εμ = εμ}} (y , t) and _{^{σμ = σμ}} (y , t) over the RVE, considering all points y of the RVE related to the point x . Thus, for an arbitrary instant t the following expressions are defined:

$\epsilon (x,t)=\frac{1}{V_{\mu }}{\displaystyle \underset{{\Omega }_{\mu }}{\int }{\epsilon }_{\mu }(y,t)dV}$ (1)

$\sigma (x,t)=\frac{1}{V_{\mu }}{\displaystyle \underset{{\Omega }_{\mu }}{\int }{\sigma }_{\mu }(y,t)dV}$ (2)

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 _^Cep , as follows:

$C^{ep}\left(x,t\right)=\frac{\partial \sigma \left(x,t\right)}{\partial \epsilon \left(x,t\right)}=\frac{\frac{1}{V_{\mu }}{\displaystyle \underset{{\Omega }_{{\mu }^S}}{\int }\partial {\sigma }_{\mu }\left(y,t\right)dV}}{\partial \epsilon \left(x,t\right)}=\frac{\frac{1}{V_{\mu }}{\displaystyle \underset{{\Omega }_{{\mu }^S}}{\int }\partial f_y\left({\epsilon }_{\mu }\left(y,t\right)\right)dV}}{\partial \epsilon \left(x,t\right)}$ (3)

Observe that after solving the RVE equilibrium problem, presented afterwards in this paper, the microscopic field of stress _^σμ is known. Then, the homogenized values σ (x, t ) and_^Cep (x, t ) can be computed from Equations (2) and (3).

Besides, the microscopic stress can be written in terms of the microscopic strain, as follows:

${\sigma }_{\mu }(y,t)=f_y({\epsilon }_{\mu }(y,t))$ (4)

where _^fy is the constitutive functional. In this work, for the triangular elements defined in the matrix (see section 3), the constitutive functional _^fy 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 _^εμ can be written in terms of the microscopic displacement filed _^uμ of the RVE, as follows:

$${\epsilon }_{\mu }(y,t)={\nabla }^{S}U{u}_{\mu }(y,t)$$ (5)

where ∇^S is the symmetric gradient operator of the displacement field u .

Without loss of generality, the microscopic displacement filed _^uu can be defined as the sum of three parts:

$u_{\mu }(y,t)=u(x,t)+{\overline{u}}_{\mu }(y,t)+{\tilde{u}}_{\mu }(y,t)$ (6)

being the first one a constant representing a rigid body motion coincident to the macroscopic displacement u (x, t ) related to the point x , the second one is obtained from the macroscopic strain ε as follows:

${\overline{u}}_{\mu }(y,t):=\epsilon (x,t)y$ (7)

which varies linearly with the coordinate y , and the last part is a displacement fluctuation field_^ũμ (y , t).

Thus, Equation (6) can be written as:

$u_{\mu }(y,t)=\epsilon (x,t)y+{\tilde{u}}_{\mu }(y,t)$ (8)

as the macroscopic displacement u (x, t ) 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)u (x, t ) has been disregarded, see (Giusti 2009).

In Equation (8) the part εy varies linearly with y resulting from the multiplication of the macroscopic strain ε of the RVE, which is constant, by the coordinates of the point y . In the case of having uniform microscopic displacement_^εμ , the displacement fluctuation _^ũμ is null. In the RVE the following relations for the microscopic strain _^εμ and the microscopic strain fluctuation ${\tilde{\epsilon }}_{\mu }$ have to be satisfied:

${\epsilon }_{\mu }={\nabla }^S{u}_{\mu }(y,t)$ (9)

${\tilde{\epsilon }}_{\mu }={\nabla }^S{\tilde{u}}_{\mu }(y,t)$ (10)

Considering Equations (8) to (10) the microscopic strain can also be written as:

${\epsilon }_{\mu }(y,t)=\epsilon (x,t)+{\tilde{\epsilon }}_{\mu }(y,t)$ (11)

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:

$${\dot{\epsilon }}_{\mu }(y,t)={\nabla }^S{\dot{u}}_{\mu }=\dot{\epsilon }(x,t)+{\dot{\tilde{\epsilon }}}_{\mu }(y,t)\forall {\dot{\tilde{u}}}_{\mu }\in {\nu }_{\mu }$$ (12)

where _^νμ 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:

${\displaystyle \underset{{\Omega }_{\mu }^S}{\int }{f}_y(\epsilon (x,t)+{\nabla }^S{\tilde{u}}_{\mu }(y,t)):{\nabla }^S\eta dV}=0\forall \eta \in {\nu }_{\mu }$ (13)

Finally the formulation is completed by the appropriated choice of the space _^νμ , 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 ε , finding the field _^ũμ ∊ _^νμ such that for each instant t, the Equation (13) is satisfied. As η is an arbitrary field, after the RVE domain discretization into finite elements, whose domain is referred as Ω_^hμ , an incremental microscopic equilibrium equation must hold for a load increment in time ∆_^tn =_{^tn +1 ^{- tn}} and a domain discretization h, finding the displacement fluctuation_{^ũμ ( n +1)} = _{^ũμ ( n )}+ ∆_{^ũμ ( n )}. If the load increment is non-linear, Equation (13) is solved by applying the Newton-Raphson Method which consists of finding the fluctuation correction ^{_δũi +1} _μ for iteration i+1, such that:

$F^i+K^i\delta {\tilde{u}}_{\mu }^{i+1}=0$ (14)

where F is the force vector and K the tangent stiffness matrix of the RVE. After computing the correction ^{_δũi +1} _μ defined in Equation (14), the next step is to obtain the displacement fluctuation field to be considered at iteration i+1 given by: ^{_ũi +1 _{μ = ũiμ + δũi} +1} _μ .

Accordingly (Pituba et al. 2016), in order to define F and K , let us denote _^Be as the element strain-displacement matrix, _^Ne as the number of triangle finite elements, _^Deμ as the constitutive tangent modulus of the element e, _^Kef as the tangent stiffness matrix of the cohesive contact finite element _^ef , _^Fintef as the internal force vector of the cohesive contact finite element and _^Nf as the number of the cohesive contact finite elements. Thus, F and K can be written as:

$F^i=\stackrel{{\Omega }_{\mu }^k}{{\displaystyle \int }}{B}^T{f}_y\left({\epsilon }_{n+1}+B{\tilde{u}}_{\mu }^i\right)dV+\stackrel{{\Omega }_{\mu }^k}{{\displaystyle \int }}{F}_f^{in{t}^{\left(i\right)}}dV={\displaystyle \sum }_{e=1}^{N_e}{B}_e^T{\sigma }_{\mu }^{e\left(i\right)}{V}_e+{\displaystyle \sum }_{e_f=1}^{N_f}{F}_{ef}^{in{t}^{\left(i\right)}}$ (15)

$$K^i=\left[\stackrel{{\Omega }_{\mu }^k}{{\displaystyle \int }}{B}^T{D}_{\mu }^{i}BdV\right]+\left[\stackrel{{\Omega }_{\mu }^k}{{\displaystyle \int }}\frac{d{F}_f^{int\left(i\right)}}{d{u}_{\mu }}dV\right]={\displaystyle \sum }_{e=1}^{N_e}{B}_e^T{D}_{\mu }^{e\left(i\right)}{B}_e{V}_e+{\displaystyle \sum }_{e_f=1}^{N_f}{K}_{ef}^{\left(i\right)}$$ (16)

$D_{\mu }^i=\left({\frac{d{f}_y}{d{\epsilon }_{\mu }}|}_{{\epsilon }_{\mu }={\epsilon }^{i+1}+B{\tilde{u}}_{\mu }^{i+1}}\right)$ (17)

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) Ω_{^{μ =}} Ω_^Sμ ⋃ Ω_^υμ , resulting into:

$\sigma =\sigma (x,t)=\frac{1}{{\text{V}}_{\mu }}{\displaystyle \underset{{\Omega }_{\mu }^S}{\int }{\sigma }_{\mu }(y,t)dV+\frac{1}{V_{\mu }}}{\displaystyle \underset{{\Omega }_{\mu }^v}{\int }{\sigma }_{\mu }(y,t)dV}$ (18)

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 whose normal direction is n ⁺ _i , must correspond to an equal side Γ^- _i with normal direction n ^- _i , being n ⁺ _i .= - n ^- _i Similarly, for each point y ⁺ defined on Γ⁺ _i must exist a point y ^- on the side Γ^- _i .To have periodic displacement fluctuation on the boundary, for every pair of points (y ⁺, y ^-) the following relation must be verified:

$${\tilde{u}}_{\mu }(y^{+},t)={\tilde{u}}_{\mu }(y^{-},t)\forall \left\{y^{+},y^{-}\right\}\in \partial {\Omega }_{\mu }$$ (19)

3 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.

Figure 2
Constitutive and FEM Modeling of the RVEs: (A) RVE with an inclusion; (B) RVE with a void.

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.

Figure 3
Contact fracture finite element developed in (Pituba et al. 2016).

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:

$$\varphi =\varphi ({\delta }_n,{\delta }_S,q)$$ (20)

where, _^δn is the normal opening displacement due to mode I; _^δS 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 =|δS|}} , therefore the behavior has an isotropic characteristic and the cohesive law is written introducing an effective opening displacement expressed by:

$$\delta =\sqrt{{\beta }^2{\delta }_S{}^2+{\delta }_n{}^2}$$ (21)

The parameter β assumes different values (from 0 to 1) establishing a weight ratio between the sliding and normal directions. The ϕ free energy potential depends of δ, and the cohesive law is expressed as:

$$t=\frac{t}{\delta }({\beta }^2{\delta }_S+{\delta }_{n}n)$$ (22)

where, n is the unit normal to the cohesive surface; _^δS 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:

$$\varphi =e{\sigma }_c{\delta }_c\left[1-e^{\left[-\left(1+\frac{\delta }{{\delta }_c}\right)\right]}\right]$$ (23)

where the law for the scalar effective traction for the loading cases is obtained from Equation (23) as:

$$t=\frac{\partial \varphi }{\partial \delta }={\sigma }_c{e}^{-\delta /{\delta }_c}if\delta ={\delta }_{\max }and\dot{\delta }\ge 0$$ (24)

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:

$$t=\frac{t_{\max }}{{\delta }_{\max }}\delta if\delta <{\delta }_{\max }or\dot{\delta }<0$$ (25)

Where e is the e-number, _^σc is the maximum tension cohesive normal traction and _^δc is a characteristic opening displacement that indicates a critical opening. Thus, β, _^σc and _^δc are parameters of the cohesive model. Besides, $\dot{\delta }$ 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:

$G_C={\displaystyle \underset{0}{\overset{R}{\int }}t.{\delta }_{\mathrm{,1}}}d{x}_1$ (26)

where R is the cohesive zone length. The Equation (26) can also be defined as:

$G_C={\displaystyle \underset{0}{\overset{\infty }{\int }}t.{\delta }_{\mathrm{,1}}}d{x}_1={\varphi }_{\infty }$ (27)

For the cohesive law presented in this work, using Equation (24), the critical energy released rate is given by:

$G_C=e{\sigma }_C{\delta }_C$ (28)

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 ). 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.

$t={\lambda }_p\delta se{\lambda }_P\delta \le {\sigma }_c$ (29)

4 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. Other analyses are performed considering porous microstructures with ductile matrix considering different void volume fractions and different void distributions. In all analyses, the periodic boundary fluctuations have been used.

Initially, the study is focused on the hydrostatic part of the homogenized stress tensor p and the homogenized effective stress of von Mises model denoted by q. Consider now, the normalized values of p’- p’ related to yielding stress_^σy , where the followings relationships are employed to obtain the numerical results:

$$p\text{'}=\frac{1}{{\sigma }_y}\left(\frac{{\sigma }_{xx}+{\sigma }_{yy}}{2}\right)$$ (30)

$$q\text{'}=\frac{1}{{\sigma }_y}\sqrt{\frac{3}{2}\left[{\left(\frac{{\sigma }_{xx}-{\sigma }_{yy}}{2}\right)}^2+{\left(\frac{{\sigma }_{yy}-{\sigma }_{xx}}{2}\right)}^2+2{\sigma }_{xy}{}^2\right]}$$ (31)

where: _^σxx represents the homogenized stress component in the direction x; _^σyy represents the homogenized stress component in the direction y; and _^σxy is the shear homogenized stress.

The numerical responses represented by curves of p’ and p’ versus normalized increments of imposed macroscopic strains have been plotted for each RVE. Soon after, the maximum values for p’ and p’ obtained in each analysis have been used for the construction of yield surfaces for the materials.

4.1 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 β = 0 (considering only normal opening of the crack surface) and β =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_{^{σy =}} 240 MPa. Also, a Young´s modulus E = 200 GPa and Poisson’s ratio ν = 0.3 have been assigned. To the inclusions, an elastic behavior has been assumed considering a Young´s modulus E = 1100 GPa and Poisson’s ratio ν = 0.25 in order to simulate a rigid inclusion. For the interface zone, the following parameters have been assigned: λ_p = 2x10¹², _^σc = 24 MPa, _^δc = 0.2mm and β = 0 or β = 0.707.

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:

$$\epsilon =\frac{\alpha }{100}\left[\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right]+\frac{\sqrt{1-{\alpha }^2}}{100}\left[\begin{array}{cc}0& \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}& 0\end{array}\right]$$ (32)

where α 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).

Figure 4
Numerical results considering the components q’ and p’ for α = 1.0.

Figure 5
Numerical results considering the components q’ and p’ for α = 0.8.

Figure 6
Numerical results considering the components q’ and p’ for α = 0.6

Figure 7
The x-direction stress distribution (MPa) over the RVE for positive values of load factor.

Figure 8
Numerical results considering the components q’ and p’ for α = 1.0.

Figure 9
Numerical results considering the components q’ and p’ for α = 0.8.

Figure 10
Numerical results considering the components q’ and p’ for α = 0.6.

Therefore, the normal strains (_^εxx and _^εyy ) and/or distortional strain (γ_xy ) vary accordingly with adopted value for α. For the analyses focusing on behavior of the normalized effective stress q’ and pressure p’, positive (α =1.0, α =0.8 and α =0.6) and negative (α =-1.0, α =-0.8 and α =-0.6) values for a are used. Besides, the analysis with α =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.

Figures 4, 5 and 6 present the results considering the components q’ and p’ for positive values of load factor.

Figure 7 presents the x-direction stress distribution over the RVE for positive values of load factor. Note that the use of a unique stress scale for the entire Figure 7 leads to some perturbations for visualization proposes of the stress distribution for α = 1.0 due to the large values of stress scale. However, the stress distribution is considered symmetric as expected for α = 1.0.

The results related to the components of q’ and p’ for compressive load factors (negative values for α) are shown in Figures 8, 9 and 10.

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 α = -1.0.

Figure 11
The x-direction stress distribution (MPa) over the RVE for negative values of load factor.

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).

$$q\text{'}=\frac{1}{1+3f+24f^6}\sqrt{1+f^2-2f\cosh \left(\sqrt{3}p\text{'}\right)}$$ (33)

where f represents the void volume fraction (in this case, f = 10%) and p’ and q’ are given by Equations (30) and (31), respectively.

Figure 12
Proposed yield surfaces considering the components q’ and p’ for four different microstructures: RVE with void, RVE with bonded inclusion, RVE with debonded inclusion (β = 0); RVE with debonded inclusion (β = 0.707).

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 β = 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 (β = 0) submitted to expansive and compressive loadings, respectively. Figures 16 and 17 are about the consideration with β = 0.707. Figure 18 presents the proposed equation for porous ductile materials.

Figure 13
Approximate function for yield surface for reinforced microstructures with perfectly bonded inclusion.

Figure 14
Approximate function for yield surface for reinforced microstructures considering phase debonding (β = 0) - expansive loadings.

Figure 15
Approximate function for yield surface for reinforced microstructures considering phase debonding (β = 0) - compressive loadings.

Figure 16
Approximate function for yield surface for reinforced microstructures considering phase debonding (β = 0.707) - expansive loadings.

Figure 17
Approximate function for yield surface for reinforced microstructures considering phase debonding (β = 0.707) - compressive loadings.

Figure 18
Approximate function for yield surface for microstructures with voids.

4.2 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 distributed 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 =}} 240 MPa. Also, a Young´s modulus E = 200 GPa and Poisson’s ratio ν = 0.3 have been assigned. It has been assumed plane strain conditions in small strain regime considering the following relation for the imposed macroscopic strain:

$$\epsilon =\frac{\alpha }{500}\left[\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right]+\frac{\sqrt{1-{\alpha }^2}}{500}\left[\begin{array}{cc}0& \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}& 0\end{array}\right]$$ (34)

where α 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).

Figure 19
Numerical results considering the components q’ and p’ for α = 1.0.

Figure 20
Numerical results considering the components q’ and p’ for α = 0.8.

Figure 21
Numerical results considering the components q’ and p’ for α = 0.6.

Figure 22
Numerical results considering the components q’ and p’ for α = 0.0.

Figure 23
The x-direction stress distribution (MPa) over the RVE for positive and null values of load factor.

Figure 24
Numerical results considering the components q’ and p’ for α = -1.0.

Figure 25
Numerical results considering the components q’ and p’ for α = -0.8.

Figure 26
Numerical results considering the components q’ and p’ for α = -0.6.

For the analyses focusing on behavior of the normalized effective stress q’ and pressure p’, positive (α = 1.0, α = 0.8 and α = 0.6) and negative (α = -1.0, α = -0.8 and α = -0.6) values for α are used as well as α = 0.0. Figures 19, 20, 21 and 22 present the results considering the components q’ and p’ for positive and null values of load factor.

Figure 23 presents the x-direction stress distribution over the RVE for positive and null values of load factor.

The results related to the components q’ and p’ for compressive load factors (negative values for α) are shown in Figures 24, 25 and 26.

Figure 27 presents the x-direction stress distribution over the RVE for negative values of load factor.

Figure 27
The x-direction stress distribution (MPa) over the RVE for negative values of load factor.

As explained in section 4.1, more analyses were performed using others values for α in order to propose the yield surfaces considering microstructures with randomly distributed voids and concentrated void presented in Figure 28.

Figure 28
Proposed yield surfaces considering the components q’ and p’for four different microstructures: RVEs with concentrated void containing 5% and 10% porosities, RVEs with randomly distributed voids containing 5% and 10% porosities.

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, 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 α = 1.0 and α = -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:

$$q\text{'}=-\frac{0.0347}{A_1}p{\text{'}}^6+\frac{0.0008}{A_2}p{\text{'}}^5+\frac{0.0777}{A_3}p{\text{'}}^4+\frac{0.0004}{A_4}p{\text{'}}^3-\frac{0.296}{A_5}p{\text{'}}^2-\frac{0.0033}{A_6}p\text{'}+\frac{0.683}{A_7}$$ (35)

Where _^Ai is a parameter to fit the curves considering centered (_{^{Ai =}} 1.0 for i = 1...7) or distributed voids

(A ₁ =1.33; A ₂ =-0.1A ₁; A ₃ =1.05A ₁; A ₄ =0.02A ₁; A ₅ =0.85A ₁; A ₆ =0.58A ₁; A ₇ =0.74A ₁).

For 10% porosity, the proposed equation is given by:

$$q\text{'}=-\frac{0.1332}{B_1}p{\text{'}}^6-\frac{0.0004}{B_2}p{\text{'}}^5+\frac{0.2104}{B_3}p{\text{'}}^4+\frac{0.0008}{B_4}p{\text{'}}^3-\frac{0.4375}{B_5}p{\text{'}}^2-\frac{0.0009}{B_6}p\text{'}+\frac{0.5794}{B_7}$$ (36)

Where _^Bi is a parameter to fit the curves considering centered (_{^{Bi =}} 1.0 for i = 1...7) or distributed voids

(B ₁ = 1.22; B ₂ = 0.01B ₁; B ₃ = 0.94B ₁; B ₄ = 0.01B ₁; B ₅ = 0.84B ₁; B ₆ = 0.12B ₁; B ₇ = 0.79B ₁).

5 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 fracture 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.

Acknowledgements

The authors wish to thank to CNPq (National Council for Scientific and Technological Development for the financial support).

References

• Abadi, M.T. (2012). Effects of microstructure on fiber-matrix debonding of Metal Matrix Composites under transverse loading. ICQNM 2012: The Sixth International Conference on Quantum, Nano and Micro Technologies, Rome, Italy.
• Azizi, R. (2012). Micromechanical modeling of damage in periodic composites using strain gradient plasticity, Engineering Fracture Mechanics 92: 101-113.
• Brünig, M., Albrecht, D., Gerke, S. (2011). Modeling of Ductile Damage and Fracture Behavior Based on Different Micromechanisms, International Journal of Damage Mechanics 20: 558-577.
• Brünig, M., Gerke, S., Hagenbrok, V. (2014). Stress-state-dependence of damage strain rate tensors caused by growth and coalescence of micro-defects, International Journal of Plasticity 63: 49-63.
• Cavancante, M. A. A., Khatam, H., Pindera, M. (2011). Homogenization of elastic-plastic periodic materials by FVDAM and FEM approaches - an assessment, Composites: Part B Engineering 42: 1713-1730.
• Chawla, N. and Shen, Y.L. (2001). Mechanical behavior of particle reinforced Metal Matrix Composites, Advanced Engineering Materials 3: 357-370.
• Chen, J.K., Huang, Z.P. e Mai, Y.-W. (2003). Constitutive relation of particulate-reinforced viscoelastic composite materials with debonded microvoids, Acta Materialia 51: 3375-3384.
• Fernandes, G.R., Pituba J.J.C., Souza Neto, E.A. (2015a). Multi-scale modelling for bending analysis of heterogenous plates by coupling BEM and FEM, Engineering Analysis with Boundary Elements 51: 1-13.
• Fernandes, G.R., Pituba J.J.C., Souza Neto, E.A. (2015b). FEM/BEM formulation for multi-scale analysis of stretched plate, Engineering Analysis with Boundary Elements 54: 47-59.
• Gărăjeu, M. and Suquet, P. (1997). Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles, Journal of the Mechanics and Physics of Solids 45: 873-902.
• Ghassemieh, E. (2002) Micro-mechanical analysis of bonding failure in a particle-filled composite, Composites Science and Technology 62: 67-82.
• Ghosh, S., Ling, Y., Majumdar, B., Kim, R. (2000). Interfacial debonding analysis in multiple fiber reinforced composites, Mechanics of Materials 32: 561-591.
• Giusti, S.M. (2009). Topological sensitivity analysis in constitutive multi-scale models, D. Sc. Thesis, LNCC/MCT, 176p.
• Giusti, S.M., Blanco, P.J.; Souza Neto, E.A., Feijóo, R.A. (2009). An assessment of the Gurson yield criterion by a computational multi-scale approach, Engineering Computations 26: 281-301.
• Gurson, A.L. (1977). Continuum theory of ductile rupture by void nucleation and growth: Part I - yield criteria and flow rules for porous ductile media, Journal of Engineering Materials and Technology, Transactions of the ASME 99: 2-15.
• Karamnejad, A., Nguyen, V.P., Sluys, L.J. (2013). A multi-scale rate dependent crack model for quasi-brittle heterogeneous materials, Engineering Fracture Mechanics 104: 96-113.
• Kattan, P. I., Voyiadjis, G. Z. (2003). A Framework for the Analysis of Damaged Composite Materials Using the Homogenization Method. Journal of the Mechanical Behavior of Materials 14: 129-148.
• Lewis, A. (2014). Making Composite Repairs to the 787, AERO 56: 04-13.
• Miehe, C. (2003). Computational micro-to-macro transitions for discretized micro-structures of heterogeneous materials at finite strains based on the minimization of averaged incremental energy, Comput. Methods Appl. Mech. Engrg.192: 559-591.
• Needleman, A., Tvergaard, V. (1987). An analysis of ductile rupture modes at a crack tip. Journal of the Mechanics and Physics of Solids 35: 151-83.
• Nguyen, V.P., Stroeven, M., Sluys, L. J. (2011). Multiscale continuous and discontinuous modeling of heterogeneous materials: a review on recent developments, Journal of Multiscale Modelling 3: 1-42.
• Papakaliatakis, G. and Karalekas, D. (2004). Study of debonding development in fibrous metal matrix composites. Proceedings of the “11th European Conference on Composite Materials (ECCM 11)”, Rhodes, Greece.
• Perić, D., Souza Neto, E.A., Feijóo, R.A., Partovi, M., Carneiro Molina, A.J. (2010). On micro-to-macro transitions for multi-scale analysis of non-linear heterogeneous materials: unified variational basis and finite element implementation, International Journal for Numerical Methods in Engineering 87: 149-170.
• Pituba, J.J.C. and Souza Neto, E.A. (2015). Modeling of unilateral effect in brittle materials by a mesoscopic scale approach, Computers and Concrete, An International Journal 15: 1-25.
• Pituba, J.J.C., Fernandes, G.R., Souza Neto, E.A. (2016). Modelling of cohesive fracture and plasticity processes in composite microstructures. Journal of Engineering Mechanics 142: 04016069.1-04016069.15.
• Reis, F.J.P. and Pires, F.M. (2013). An adaptive sub-incremental strategy for the solution of homogenization-based multi-scale problems, Comput. Methods Appl. Mech. Engrg 257: 164-182.
• Santos, W.F., Fernandes, G.R., Pituba, J.J.C. (2016). Análise da influência dos processos de plasticidade e fratura no comportamento mecânico de microestruturas de Compósitos de Matriz Metálica, Matéria (UFRJ) 21: 577-598.
• Segurado, J. and Llorca, J. (2005). A computacional micromechanics study of the effect of interface decohesion on the mechanical bahavior of composites, Acta Materialia 53: 4931-4942.
• Somer, D.D., Perić, D., Souza Neto, E.A., Dettmer, W.G. (2015). Yield surfaces of heterogeneous media with debonded inclusions, Engineering Computations 32: 1802-1813.
• Sun, L.Z., Ju, J.W., Liu, H. T. (2003). Elastoplastic of metal matrix composites with evolutionary particle deboing, Mechanic of Materials 35: 559-569.
• Tian, R., Chan, S., Tang, S., Kopacz, A.M., Wang, J.-S., Jou, H.-J., Siad, L., Lindgren, L.-E. Olson, G.B., Liu, W.K. (2010). A multiresolution continuum simulation of the ductile fracture process, Journal of the Mechanics and Physics of Solids 58: 1681-1700.
• Tvergaard, V., Needleman, A. (2006). Three dimensional microstructural effects on plane strain ductile crack growth. International Journal of Solids and Structures 43: 6165-79.
• Tvergaard, V. (1981). Influence of voids on shear band instabilities under plane strain conditions, International Journal of Fracture 17: 389-407.
• Voyiadjis, G. Z., Kattan, P. I. (1993). Local Approach to Damage in ElastoPlastic Metal Matrix Composites. International Journal of Damage Mechanics 2: 92-114.
• Zhang, M.H. and Chen, J.K. (2012). Analysis of interfacial fracture strength of an inclusion in a polymeric composite considering cohesive force, Computational Materials Science 61: 6-11.

Publication Dates

• Publication in this collection
  16 Aug 2017

History

• Received
  27 Feb 2017
• Reviewed
  07 May 2017
• Accepted
  23 May 2017