On void shape effects of periodic elasto-plastic materials subjected to cyclic loading

This work investigates the effect of cyclic loading upon periodic elasto-plastic porous materials. The aim is to understand the evolution of the underlying microstructure, described here with a single void embedded in a cubic unit-cell. Periodic ﬁ nite element (FEM) calculations are carried out under a ﬁ nite strain deformation process keeping the absolute value of the stress triaxiality and the Lode angle con-stant during the cycle. As a result of the applied loading conditions, the void geometry, both volume and shape, change signi ﬁ cantly leading to porosity and void shape ratcheting. The void shape becomes non-spherical from the very ﬁ rst cycle leading to a markedly asymmetric cyclic response of the material. This, in turn, results in an observed maximum stress as a function of the number of cycles. In addition, even though the average applied strains are relatively small, the local strains near the void surface increase signi ﬁ cantly inducing a signi ﬁ cant localization of the deformation. Finally, several initial void shape con ﬁ gurations are also considered. In the majority of the cases studied, the void evolves into a crack-type shape in the direction of the minimum absolute stress. This, in turn, is consistent with a con ﬁ guration corresponding to a crack subjected to a mode I cyclic loading.


Introduction
Cyclic loading of metallic materials has always attracted a lot of attention in the scientific and industrial community due to its impact in the numerous low and high cycle fatigue applications. A priori, the response of a metal under cyclic loading conditions is a multi-scale problem and several mechanisms at the microscale (e.g., dislocation structures at the scale of 1e10 mm) and meso-scale (e.g., hard precipitates and pores or cracks) lead to damage at different scales and finally initiate fatigue of the material at the macroscale. This makes the analysis and modeling of cyclic loading a very tedious work in the sense that to-date it is very difficult to propose a micromechanics-based model that spans the fatigue mechanisms in all the scales. Nonetheless, a lot of studies have been made at several length scales and many of the mechanisms in cyclic loadings have been identified.
In the high cycle regime that plasticity is strongly confined, Dang Van (1971) (see also Papadopoulos (1987) and Constantinescu et al. (2003)) proposed pressure dependent fatigue criteria based on a combined homogenization and phenomenological approach. Numerical modeling of the two scale approach (Bertolino et al., 2007;Hofmann et al., 2009;Guerchais et al., 2014) shows that the pressure is the footprint of the residual stress created by the localized plasticity at the grains scale. In turn, in the low cycle regime, where plasticity is more spread at numerous grains, only phenomenological and statistically-based fatigue criteria have been proposed (see for example Amiable et al. (2006a,b) and Tabibian et al. (2013b)). Mean stress effects on high strength steels have been discussed already by Koh and Stephens (1991) and Kondo et al. (2003). Recently, Morel and Bastard (2003) and Maitournam et al. (2011) for high-cycle fatigue and Tabibian et al. (2013a) for low-cycle fatigue have experimentally shown the effect of multiaxial loading effects, and particularly of pressure, upon the cyclic response of steels and aluminum alloys and have introduced similar parameters for taking into account the mean stress effect. The fact that fatigue criteria, which are based on strong experimental evidence, include pressure dependence, is in striking contradiction with the common modeling of such materials using pressure-independent plasticity laws, both crystal plasticity and isotropic plasticity with or without kinematic hardening.
More specifically, advances in imaging techniques (SEM and tomography) have revealed the presence of voids in metals. In a recent study, Limodin et al. (2014) and Wang et al. (2014) have obtained 3D tomographic images (see Fig. 1) for aluminum alloys obtained by lost-foam-casting fabrication techniques. As shown in Fig. 1c, voids of several sizes and shapes (spherical, ellipsoidal but also non-canonical) are observed. In the same figure, the different colors indicate families of voids of similar shapes (but not size or orientation), while two representative voids are pointed out by arrows. One is an almost spherical void (with aspect ratios w 1 z w 2 z 1 and the other is a prolate void with average shape w 1 z w 2 z 2 (see Fig. 5a for a detailed definition of void aspect ratios). It is further noted that the size of the voids ranges between 50À500 mm with grain size in the order of 50À100 mm.
Motivated by the tomographic image in Fig. 1, a more physicsbased way to include such pressure dependence at the material level is the use of void microstructures 1 embedded in an otherwise plastically incompressible matrix phase. These voids or cracks, which could be present in the material ab initio or be nucleated around precipitates and particles in the course of deformation (Essmann et al., 1981), could be smaller, equal or larger than the size of the grains. The voids therefore can have sizes from a few microns (e.g., 1e10 mm) to hundreds of microns (larger than 200 mm). The presence of the pore, in turn, gives rise to a compressible response of the composite material in the plastic region, and hence to pressure dependence. This is achieved by transforming the mean stress applied at the macroscale to a local shearing of the matrix material near the void surface.
As a consequence, the material is now viewed as a two-phase composite system comprising the void phase and the matrix phase. The matrix phase, depending on which scale we refer to, could be the grain or an ensemble of grains. In the first case, that the voids lie inside the grain or the cracks are in the scale of grains, it was shown that the strain-gradient effects with (Niordson and Legarth (2010), Vernerey et al. (2007) or without (Monchiet and Kondo (2013)) crystal plasticity can have a significant impact on material response Miehe et al. (1999); Watanabe et al. (2010). In the second case, where the voids are larger than the size of the grains, standard isotropic plasticity could be used to model the response of the material, which is the case in the present study.
While the idea of using porous materials to study ductile fracture of metallic materials subjected to monotonic loading conditions has been used extensively in the literature, very few studies have been carried out in the domain of cyclic response of such materials.
Specifically, in the context of monotonic loading conditions, nonlinear homogenization models and micromechanical models (such as the well-known model of Gurson (1977) for elasto-plastic porous materials) have been used for the prediction of pressure dependent material behavior. In the context of nonlinear homogenization techniques, Ponte Castañeda (1991) (see also Michel and Suquet (1992) for a parallel development using a different approach) has proposed a linear comparison composite method initially applied to isotropic porous materials. In a later stage, these techniques have been extended to include in an accurate manner general ellipsoidal void shapes (see for instance the recent works of Danas and Ponte Castañeda (2009a) and Danas and Aravas (2012)). In a parallel development, extensions of the Gurson model either for isotropic void shapes (see for instance Tvergaard (1984) and Mear and Hutchinson (1985)) or spheroidal (Gologanu et al., 1993;Benzerga and Besson, 2001;Monchiet et al., 2006) and ellipsoidal (Madou and Leblond, 2013) void shapes have been proposed. Such material systems have also been analyzed very early using numerical finite element methods (see for instance the seminal work of Koplik and Needleman (1988)), and are still addressed by recent works (see for instance Tvergaard (2011);Tvergaard. (2011, 2012); Tekoglu et al. (2012)). In these works, the pressure dependence has been studied via the stress triaxiality, denoted here by X S and defined as the ratio of the mean stress to the von Mises equivalent or effective deviatoric stress. More recently, the Lode angle q, which is directly related to the third invariant of the deviatoric stress tensor, has been identified experimentally (Barsoum and Faleskog, 2007a) as an important loading parameter, especially at low stress triaxialities. Numerical simulations (Zhang et al., 2001) and analytical micromechanical models have been proposed in this regard (Nahshon and Hutchinson, 2008;Danas and Ponte Castañeda, 2012).
Nevertheless, most of the above studies have been carried out in the context of monotonic loading conditions. Even though the material is initially the same, the evolution of the void size and shape in cyclic loading conditions is markedly different than in the context of monotonic loadings. Yet, much less has been done in the context of cyclic loading conditions. Specifically, Monchiet et al. (2008b) have used a micromechanical model for porous materials to explain the mean stress effect in high cycle fatigue. Furthermore, Devaux et al. (1997), Besson and Guillemer-Neel (2003) and Rabold and Kuna (2005) have explored numerically the cyclic response of porous materials at small and moderate number of cycles with a main emphasis on axisymmetric loading states at large strains. Their analysis has mainly focused on the prediction of porosity ratcheting, whereby the underlying void shape changes have not been studied in detail. In a similar study, Ristinmaa (1997) has (c) Voids of sizes ranging between 50À500 mm. The two arrows show voids with average aspect ratios w 1 z w 2 z 2 (prolate ellipsoidal shape) and w 1 z w 2 z 1 (average spherical shape). Courtesy of N. Limodin and E. Charkaluk. 1 To avoid any misunderstanding with the different communities in mechanics and material science, we precise here that the word "microstructure" refers, henceforth, to the voids. carried out finite-strain unit-cell computations, but for a small number of cycles, concluding that void shape effects have very little effect on porosity ratcheting.
In a slightly different context, Pirondi et al. (2006) and Hommel and Meschke (2010) have used the Gurson (1977) and the Leblond et al. (1995) models, respectively, to investigate the low cycle fatigue response of metallic structures. Rather interestingly, the later found that by including void shape effects (contrary to the first who used Gurson model that includes no void shape effects) could dramatically improve their predictions, even though a large number of additional fitting parameters had to be used. Similar observations regarding the importance of void shape effects upon the cyclic response of porous materials have also been made recently by Carpiuc (2012) who used the model of Danas and Aravas (2012), which includes general ellipsoidal void shapes and orientations. In that study, it was found that void shape changes tend to accumulate in each cycle thus leading to ellipsoidal void shapes and consequently to porosity ratcheting, contrary to the Gurson model that predicts no porosity ratcheting (see Devaux et al. (1997) for more details). Nonetheless, all these homogenization and micromechanical models contain only partial information about the void shape changes (up to a perfect ellipsoidal shape) and as we will see in the following they need to be re-assessed first via numerical calculations such as the   3. (a) Schematic representation of the principal stress (s 1 ,s 2 ,s 3 ) cartesian system and the (s m ,s eq ,q) cylindrical system. (b) Components of the normalized stress 3s i /2s eq , i ¼ 1,2,3 as a function of the Lode angle q in the case of X S ¼ 3 Fig. 4. Schematic explanation of the application of the cyclic loading and the corresponding qualitative values of (a) the applied displacement u, (b) the applied stress triaxiality X S and (c) the applied Lode angle q as a function of time for one cycle.
present study and ultimately be used to compare with experiments.

Scope of this study
The scope of the present study is to investigate the effect of cyclic loading conditions and finite deformations upon microstructure evolution and material softening/hardening using finite element (FEM) periodic unit-cell calculations with 3D geometry at small and large number of cycles. The matrix material is described by isotropic J 2 plasticity considering the case, described previously, of the voids being much larger than the grain size but smaller than the specimen size. Nonetheless, the results obtained in the present study could still be valid, at least in a qualitative manner, in the case of materials with large, but not strongly anisotropic (e.g., FCC) grains comprising intragranular voids that are not of a nanometer size (Monchiet and Kondo, 2013).
Furthermore, it is worth to impress upon the fact that in the fatigue community, a large majority of studies in the context of cyclic loadings is done using small strain calculations, in order to increase numerical efficiency, and therefore neglect any changes of the underlying microstructure (including void shape effects). As we will see in this work, however, the local strains can be in excess of 50% due to strong localization of the deformation around voids, even if the overall applied strains are small (in the order of a few percent). For that reason and in agreement with the aforementioned micromechanical studies (Devaux et al., 1997;Besson and Guillemer-Neel, 2003;Rabold and Kuna, 2005), it is also critical that a finite deformation analysis is carried out in the present work.
Specifically, in Section 2, we describe the unit-cell geometry and the applied loading states at the scale of the microstructure, defined here by the presence of a void in an homogeneous elastoplastic matrix. Furthermore, we identify the variables used to analyze the cyclic response of the periodic porous medium and the void geometry changes. Next, in Section 3, we present the cyclic response of the unit-cell at small and large number of cycles where we identify the principal micro-deformation mechanisms that lead to an overall softening of the porous material. In the following, in Section 4, we carry out a parametric study in order to investigate the effect of the loading and the initial void shape, respectively, upon the cyclic response of the unit-cell. It should be mentioned here that the above described sections are devoted to a plastically incompressible matrix phase with purely isotropic hardening. However, in Section 5 preliminary calculations with coupled nonlinear isotropic-kinematic hardening will also be considered showing similar qualitative characteristics with the purely isotropic hardening case. Finally, we conclude with a brief discussion of the main results and perspectives of this study.

Problem formulation
In this section, we define a periodic porous medium with cubic unit-cell geometry as well as the loading conditions used in this study. The interest in this work is the analysis of a cubic periodic unit-cell comprising a spherical void positioned at the center. The unit-cell is subjected to periodic average cyclic loading conditions with a constant amplitude of average stress triaxiality and average Lode angle. A critical aspect of the present study is the use of a finite strain analysis contrary to the more common small strain analysis used when studying cyclic loading conditions. This will allow for the evolution of the void geometry due to the local large strains at the current configuration. In this regard, one needs to identify the relevant parameters that are necessary to describe the evolution of the void geometry, in general.

Geometry of the unit-cell
In order to set the stage of the following analysis we first attempt to separate the relevant length-scales of the problem at hand. In this regard, let us consider a three-dimensional specimen as shown in Fig. 2a. It is common practice to consider that the material is homogeneous and is described by phenomenological constitutive laws (e.g., J 2 plasticity, anisotropic plasticity, etc). In reality, however, this specimen is heterogenous and in several cases of practical interest comprises defects, e.g., impurities, cracks and/ or pores at the micron scale (see Fig. 1) at a scale which is bigger than the scale of the size of the grains but much smaller than the size of the specimen and the scale of variation of the externally applied loading conditions. Furthermore, such a specimen is usually subjected to general loading conditions (such as traction, torsion or a combination of both, etc), whereas at the local level, one finds a rather complicated stress and strain loading state due to the nonlinear macroscopic geometry, which involves both shears and hydrostatic stress states. As already stated in the previous section, it has long been acknowledged that the effect of pressure is primordial in cyclic loading conditions, yet such materials are modeled very often as plastically incompressible. In this work, in order to include (physically) this pressure dependence, we consider that the underlying material that makes up the notched specimen is a periodic medium comprising initially spherical voids distributed with cubic symmetry as shown in Fig. 2b, whereas the two scales e specimen scale and material scale e are well separated. Even though the choice of a cubic symmetry is an idealized choice, the use of a small initial void volume fraction (i.e., <5%) together with overall low strains (i.e., 1%) allows for a sufficiently general qualitative analysis without any significant interactions between neighboring voids, as we will see in the following, at least until very late where void shapes evolve significantly and mesh distortion is prohibitive for further numerical analysis. In addition, the use of a periodic medium allows for a full-field numerical analysis of the material response due to the fact that only one cubic unit-cell ( Fig. 2c) with appropriately defined periodic boundary conditions is needed (Michel et al., 1999). This single unit-cell can then generate by periodic repetition the entire microstructure of the composite (Fig. 2b).
More specifically, the presence of voids will immediately give rise to an average hydrostatic stress dependence of the periodic medium at the plastic range since the voids are compressible and hence the average (plastic) strain in the unit-cell will have a nonzero hydrostatic component (i.e., it exhibits compressible plasticity). In order to keep the analysis tractable as well as to simplify our numerical calculations, we will further restrict attention to only triaxial loading states aligned with the symmetries of the unit-cell, thus analyzing only one-eighth of the cube, as shown in Fig. 2d. This implies that the underlying void will evolve in volume and in shape when finite strains are applied at the level of the unit-cell, but not in orientation. 2 Finally, in this study, we attempt to make no direct coupling between the several scales, i.e., the specimen scale (Fig. 2a) and the material scale ( Fig. 2b) but mainly to understand the effect of a triaxial stress state upon the cyclic response of a periodic porous material. This, of course, implies further that we cannot carry out any direct comparison with experiments since these require the use of a geometry, such as the one in Fig. 2a.
On the other hand, this analysis has as a focus to analyze and understand the basic microstructural deformation mechanisms so that (less time-consuming) analytical homogenization models for porous materials (Monchiet et al., 2008a;Danas and Aravas, 2012;Madou and Leblond, 2013) and phenomenological pressuredependent criteria (Dang Van, 1971;Constantinescu et al., 2003;Monchiet et al., 2008b) can be re-assessed.

Periodic boundary conditions and cyclic loads
We consider a cubic unit-cell occupying a volume V with side length of 2L and boundary vV that comprises a spherical void at the center, as shown in Fig. 2c. The matrix phase is described by an isotropic elasto-plastic constitutive law as described later in this section.
The unit-cell is subjected to triaxial periodic boundary conditions, so that the velocity field _ u (the superposed dot denotes time derivative) can be split into an affine part D,x and a correction _ u * , i.e., (Michel et al., 1999;Miehe et al., 1999) _ The second-order tensor D, characterizing the affine part, corresponds to the average strain-rate field in the periodic medium (i.e., the actual strain-rate field of the unit-cell if it were homogeneous) and is defined via the local fieldD which admits the following decomposition.
The average strain ε at a given time t in the unit-cell is expressed by.
In turn, the average stress s (i.e., the actual stress field of the unit-cell if it were homogeneous) is defined formally in terms of the local stress fields via.
wheres satisfies the following equilibrium equations and periodic boundary conditions, i.e., divs ¼ 0 in V;s$n À # on vV: In this expression, n denotes the normal to the exterior faces of the unit cell and À# is used to denote that the traction is opposite on opposite sides of the unit-cell. Henceforth, the use of the quantities s, D and ε (and of any other quantity resulting from those) refers unambiguously to the average (or macroscopic) stress, strain-rate and strain fields in the unit-cell.
As a consequence of the presence of the void, the average response of the unit-cell depends upon the average hydrostatic pressure (or mean average stress) as well as the deviatoric part of the average stress tensor. Thus, it is useful to define at this point the average stress triaxiality, X S , and average Lode angle, q, in the unitcell as measures of the average stress state in the unit-cell, such that.
with s'¼s À s m I denoting the stress deviator and I the secondorder identity tensor. Using the definitions in equation (6), one can write the principal components of the stress field as a function of X S and q, via 3 2s eq fs 1 ;s 2 ;s 3 g¼ n cosq;Àcos The graphical illustration of the above relations is shown in Fig. 3. The (s m ,s eq ,q) coordinates define a cylindrical coordinate system oriented along the hydrostatic axis s m , as shown in Fig. 3a.
In Fig. 3b the normalized stress components, 3s i /2s eq are depicted as a function of the Lode angle q for a stress triaxiality X S ¼ 3. Note that due to the p/3 periodicity of the functions used in relation (7), the three principal stresses exhibit similar periodicity. For later use, we also show explicitly in Table 1 the order of the stress components for three representative Lode angles that will be used in the following sections. Note that in the case of an initially non-spherical shape a larger range of Lode angles (i.e., q > 60 o ) could be considered. However, for the sake of brevity, we will restrict attention only to the aforementioned Lode angles.
In the following, we focus on purely triaxial loading conditions aligned with the symmetry planes of the cubic unit-cell and the underlying void microstructure. This allows for numerical simplicity while keeping essential features of the void geometry changes. This leads to the following Dirichlet-type boundary conditions.

<
: with u * ¼ 0 on vV. This implies that the external faces of the unitcell remain straight (Michel and Suquet, 1994;Garajeu et al., 2000) and hence only 1/8 of the unit cell may be considered, as shown in Fig. 2d. In turn, the void surface is traction free.
Before we define the cyclic loading in terms of the displacements U i (t) (i ¼ 1,2,3), we discuss, first, the methodology for the Table 1 Ordering of the principal stresses for different Lode angles. The relevant directions are shown in Fig. 2d while their relevant magnitude and graphical representation is shown in Fig. 3b.
Orientation effects can readily be included by considering the entire unit-cell but this is left for a future study since this would introduce a large number of additional set of parameters to be investigated and would make the present work too lengthy.
application of a constant stress triaxiality and Lode angle during a given deformation process. This methodology has been originally proposed by Barsoum and Faleskog (2007b) (see also Dunand and Mohr (2014)) and is further discussed here for completeness of the present manuscript. As a consequence of the applied periodic boundary conditions and the symmetry of the problem at hand, the average deformation in the unit-cell is entirely described by the three displacements on its exterior surfaces (c.f. equation (8)), denoted compactly as.
U ¼ fU 1 ðtÞ; U 2 ðtÞ; U 3 ðtÞg: Recalling that the average strain-rate and stress tensors involve only three non-zero components due to the applied triaxial loading, they can be expressed in vectorial form (i.e., using the Voigt notation) as.
To proceed further, we rewrite the strain-rate tensor as.
where Q is diagonal matrix of dimension three. We, next, define an external fictitious node, 3 whose generalized force, P G , and generalized displacement, p G , vectors, respectively, take the form (see also Lin et al. (2006) for equivalent formulations) The stress state in the unit-cell is then controlled via a timedependent kinematic constraint (Michel et al., 1999) obtained by equilibrating the rate of work in the unit-cell with the rate of work done by the fictitious node on the unit-cell at time t, such that Next, in order to control the loading path in the stress space, we couple the average stress s in the unit-cell with the generalized force vector associated with the fictitious node P G via the constraint equation.
where C is a non-dimensional proper orthogonal matrix since c i (i ¼ 1,2,3) are three dimensional vectors that form an orthogonal basis set. The vectors c i (i ¼ 1,2,3) depend on the three components of the average stress s, such that The above expressions for the vectors c i (i ¼ 1,2,3) together with the definitions (7) further imply that the matrix C in equation (14) is only a function of the stress triaxiality X S and the Lode angle q but not of the equivalent stress s eq . By substitution of equations (11) and (14) in (13), one gets.
The above expression provides the kinematic constraints between the degrees of freedom corresponding to the sides of the unit-cell (i.e., U) and the degrees of freedom of the fictitious node (i.e., p G ). These nonlinear constraints are applied in the finite element software ABAQUS (ABAQUS, 2009) by use of the multipoint constraint user subroutine (MPC).
Using the above definitions, we divide each cycle in four steps. In each step of each cycle, as shown in Fig. 4a, we control the average strain in the unit-cell, by setting uðtÞ=L≡p G 1 ðtÞ. 4 This quantity, u/L, initially increases from O to A (step 1), unloads from A to B (step 2), reversely loads from B to C (step 3) and unloads from C to D (step 4) defining thus an entire cycle. Then, the average strain-rate D is evaluated such that the stress triaxiality X S and Lode angle q remain constant in each step as discussed previously and as shown in Fig. 4b and c. Note in these figures that in order to obtain full stress reversibility during the cycle, X S has to change sign and q has to jump to q þ p between A-D. For convenience, hereafter, the notations X S and q are used to denote unambiguously the absolute value of the stress triaxiality (i.e., X S ≡jX S j) and the minimum value of the Lode angle (i.e., q≡cos À1 cosq ) in each cycle AeD.
In the following calculations, the matrix phase is described by an elasto-plastic constitutive relation. The elastic part is defined via the Young's modulus E and the Poisson ratio n. In turn, standard J 2 plasticity theory is used to describe the plastic behavior of the matrix together with an isotropic strain hardening law (except in Section 5 where a nonlinear kinematic hardening law is also added) given by a power-law form, which reads Here, s 0 and ε 0 denote the initial yield stress and yield strain of the matrix material, N is the hardening exponent and ε p is the accumulated plastic strain in the matrix phase defined in the usual way. In this study, we focus on realistic values of the elastic moduli, e.g., Young's modulus, E~1000s 0 (for instance s 0~2 00 MPa corresponding to steel), Poisson's ratio, n ¼ 0.3 and hardening exponent, N ¼ 10. Nonetheless, one should point out that particularly in cyclic loading conditions, the effect of elasticity and hardening could be important as already discussed in Devaux et al. (1997), and such an analysis is detailed in Appendix A.

Evolution of void geometry
In this section, we introduce the variables used to characterize the evolution of the change in volume and shape of the void. More specifically, the porosity (i.e., the volume fraction of the void in the unit-cell) is defined as.
where V v , V m and V are the volume of the void, the matrix and the total volume of the unit-cell, respectively. Here V m is calculated as the sum of each volume element, while the unit-cell volume V is evaluated using the coordinates of the corner nodes of the cubic unit-cell since due to symmetry of the void and the purely triaxial loading conditions the external faces of the cell remain straight. As a consequence of the finite deformations considered in this study, significant changes in the pore shape are also observed. Therefore, appropriate geometrical quantities need to be introduced in order to evaluate such pore shape changes. As a first-order measure, the void shape is characterized by two aspect ratios 3) denote the current lengths of the axes of the void that intersect with the coordinate axes x i (with i ¼ 1,2,3), respectively, as shown on Fig. 5a. However, such a measure will be shown in the following to be insufficient since the void obtains markedly non-ellipsoidal shapes due to the cyclic loadings as opposed to the purely monotonic loadings where almost ellipsoidal shapes are observed Needleman, 2012, 2013). In this regard then, as an additional measure of the pore geometry change, we have also defined the ellipsoidicity ratio. This ratio has been introduced as a measure of the divergence of the void geometry from an equivalent perfect ellipsoid, as depicted in Fig. 5b. While a large number of options can be used to identify this difference, use is made here of a simple measure. First, we set the axes of the ideal ellipsoid equal to the length of the actual void axes. Then, the volume of the ideal ellipsoid, V e , will in general be different from that of the actual void V v due to the nonlinearity of the matrix phase, the interactions of the neighboring voids of the periodic composite and more importantly due to the cyclic loading conditions. Therefore, the ellipsoidicity ratio, ε l defined via gives the difference of the actual void shape from that of a perfect ellipsoid. Consequently, when the ellipsoidicity ratio takes values close to unity, the void shape remains almost an ellipsoid. The above microstructural variables will be used in the following to analyze the micromechanisms that lead to the material softening/hardening due to the applied cyclic loading conditions. It should be noted here that due to the cubic symmetry of the unitcell and the purely triaxial loading conditions no void rotations are obtained.

Cyclic response and microstructure evolution
In this section, we discuss the results obtained by the previously described loading conditions. The cyclic loading conditions are parametrized by the use of two different values of the stress triaxiality X S ¼ 2/3,3 and three different values of the Lode angle q ¼ 0 ,30 ,60 are used. 5 For the low triaxiality X S ¼ 2/3, we set the average strain amplitude u/L ¼ 5% and for the high triaxiality X S ¼ 3, we set the average strain amplitude u/L ¼ 1%. The difference in amplitudes has been introduced for convenience with the calculation time needed to observe significant void geometry changes and/or localization of the strain at certain region of the unit-cell. Two sets of computations have been carried out: (i) for a small number of cycles, e.g., 5 cycles and (ii) for a large number of cycles, e.g., z50 cycles.
Moreover, for convenience with the meshing, we use an initial porosity f ¼ 1%, which corresponds to a void radius a/L ¼ 0.2673, and a maximum of 32 Â 10 4 degrees of freedom, which leads to an average of 2.5 h computational time per cycle on a 12-cpu parallel computation. A detailed discussion of dependency of the results upon the mesh size is carried out in the Appendix B.
In order to clarify further the results in the following sections, we include Fig. 6, where for a given variable A (e.g., porosity, aspect ratios, ellipsoidicity, etc), an average quantity per cycle is defined as the arithmetic mean of the maximum and the minimum value per cycle. In the same figure, the flat, horizontal part of the curves corresponds to elastic unloading. On the other hand, the average von Mises stress per cycle is evaluated at the end of the first step of each cycle, as those are defined in the context of Fig. 4.

Small number of cycles
In Fig. 7, we consider stress-strain results for 5 cycles, a low stress triaxiality X S ¼ 2/3 and Lode angle q ¼ 0 (with u/L ¼ 5%).
Specifically, we show two normalized average principal stresses, (a) s 1 /s 0 and (b) s 2 /s 0 , in the unit-cell as functions of the average components of the strains, ε 1 and ε 2 , respectively. We focus on these two stress components since Lode angle q ¼ 0 corresponds to an axisymmetric stress state case (see Table 1) and hence the third stress component is equal to the second one. In Fig. 7a, we observe a common stress-strain cyclic response where in the first cycle a hardening is obtained in both tension and compression. This hardening tends to saturate rather fast due to the low hardening exponent used in this case, as observed by the use of the notation C1, C2, C3, C4, C5 which serve to identify each of the 5 cycles. In Fig. 7b, the second average stress-strain response is similar to the first one but with a different sign of the strain due to the plastic incompressibility of the matrix. Note however that the average plastic response of the unit-cell is not incompressible due to the presence of the void, however, due to the low stress triaxiality in this example, this is only slightly visible by noting the small asymmetry of the curve in Fig. 7b with respect to the average strain ε 2 .
Similarly, we show two normalized average principal stresses, (a) s 1 /s 0 and (b) s 2 /s 0 , in the unit-cell as functions of the average components of the strains, ε 1 and ε 2 , respectively. As observed in Fig. 8a, the material significantly softens during the cycle when ε 1 is (positive) tensile while it hardens when ε 1 is (negative) compressive. As will be seen next, this is due to the evolution of the porosity (i.e., void volume fraction), which increases for ε 1 > 0 and decreases for ε 1 < 0. This leads, in turn, to a significant asymmetry of the stress-strain response between positive and negative stress triaxialities (i.e., between tension and compression) which is markedly different than the corresponding stress-strain response at the lower stress triaxiality of X S ¼ 2/3 of the previous figure. Note that the observed asymmetry is not due to a Bauschinger effect but is strongly related to the evolution of the void geometry as will be seen in the following. This asymmetry has also been identified in the context of uniaxial yielding by Cazacu et al. (2014) as a possible cause of swift effect Swift (1947) i.e. the occurrence of inelastic axial effects. On the other hand, in Fig. 8b, we observe a strong asymmetry in the strain axis. This is due to the high triaxiality loading used in this case (X S ¼ 3), resulting to a highly compressible plastic response of the unit-cell. In Fig. 9, we discuss the evolution of the porosity f for the two afore-mentioned stress triaxialities, i.e., for (a) X S ¼ 2/3 and (b) X S ¼ 3 with q ¼ 0 as a function of the average strain ε 1 . As shown in Fig. 9a, only a minor porosity ratcheting is observed if one observes the extremities of the cyclic curves at positive strains. In addition, as already stated previously, f increases for ε 1 > 0 and decreases for ε 1 < 0, as intuitively expected. On the contrary, porosity evolution is much more significant for X S ¼ 3, as observed in Fig. 9b, where porosity grows by almost 50% (note the crossing of the f curves at ε 1 ¼ 0) after only 5 cycles. This can explain the observed asymmetry of the stress-strain curve with respect to the strain ε 2 in Fig. 8b.
At this point, it is worth noting that the observed effect of stress triaxiality upon the above-discussed stress-strain responses and 5 Recall that the notation X S ¼ 2/3 and q ¼ 0 , for instance, corresponds to jXSj ¼ 2=3 and q ¼ cos À1 cos q according to the discussion made in the context of Fig. 4. porosity ratcheting, for a matrix with purely isotropic hardening, is in full qualitative agreement with the results presented by Rabold and Kuna (2005) in the context of combined isotropic and kinematic hardening. This, further, implies that the hardening characteristics of the matrix phase affect only quantitatively, but not qualitatively, the evolution of the void volume and shape. This observation is further confirmed in Section 5. Fig. 10 shows the evolution of the void shape via the evolution of the aspect ratio w 1 for the two stress triaxialities (a) X S ¼ 2/3 and (b) X S ¼ 3 with q ¼ 0 as a function of the average strain ε 1 . Note that the second aspect ratio w 2 ¼ 1, since for q ¼ 0 the loading is axisymmetric along x 1 (i.e., js 2 j ¼ js 3 j). More specifically, in Fig. 10a, we observe that for the lower stress triaxiality X S ¼ 2/3, the void aspect ratio exhibits only a small ratcheting. Interestingly, as the number of cycles increases, the void shape tends to become oblate (i.e., w 1 > 1, a 1 < a 3 ¼ a 2 ) even though js 1 j > js 2 j ¼ js 3 j and the aspect ratio decreases in the first step of C1, as intuitively expected. However, the very small porosity and aspect ratio ratcheting in this low triaxiality case leads to an almost symmetric stress-strain curve in Fig. 7. On the other hand, in Fig. 10b, a significantly asymmetric void shape change with respect to ε 1 is observed. Moreover, as a result of the high stress triaxiality, the void tends to become oblate (w 1 > 1, a 1 < a 3 ¼ a 2 ) from the very first loading step even though js 1 j > js 2 j ¼ js 3 j. This counterintuitive effect has been attributed to the strong nonlinearity of the matrix and the high triaxiality loading (see Budiansky and Hutchinson (1980) and Fleck and  stress s 2 /s 0 as a function of the average strain ε 2 at 5 cycles in the case of u/L ¼ 5%, X S ¼ 2/3, q ¼ 0 . The notation C1, C2, C3, C4, C5 represents the first, second, third, fourth and fifth cycle, respectively. Fig. 8. (a) Normalized average stress s 1 /s 0 as a function of the average strain ε 1 and (b) normalized average stress s 2 /s 0 as a function of the average strain ε 2 at 5 cycles in the case of u/L ¼ 1%, X S ¼ 3, q ¼ 0 . The notation C1, C2, C3, C4, C5 represents the first, second, third, fourth and fifth cycle, respectively.
Hutchinson (1986)). More importantly, the oblate shape persists as the number of cycles increases. This is due to the asymmetric response of the unit-cell between ε 1 > 0 and ε 1 < 0. More specifically, we find that the rate of the void shape change is more significant in the ε 1 < 0 (i.e., in compressive loads) rather than in the ε 1 > 0 (i.e., in tensile loads). This asymmetry can be qualitatively observed by the difference in the slopes of the w 1 -ε 1 curves in the ε 1 > 0 and ε 1 < 0 regimes. In turn, this asymmetry in the void shape evolution leads to a permanent irreversible void shape change from the very first cycle (i.e., after the first cycle the void is not spherical), which in turns produces the porosity ratcheting and the asymmetric average stress-strain response in Fig. 8. At this point it is worth mentioning that the void shape irreversibility and consequently porosity ratcheting are present also for different hardening exponents (including N / ∞ which corresponds to an ideallyplastic response) not shown here for the sake of a reasonable set of parameters investigated.

Large number of cycles
In the previous section, it has been shown that the cyclic behavior of the material, and especially the porosity ratcheting as well as the asymmetric void shape evolution have a strong effect upon the average stress-strain response of the unit-cell. Thus, in order to isolate the effect of porosity evolution and the influence of the void shape change in the cyclic response of the material, a unitcell geometry with a spherically constrained void shape is further added in our study. It is further clarified here that the spherically constraint void does allow the evolution of the void volume (i.e., void growth) and porosity change but constraints the change of the void shape to remain spherical. This is achieved by imposing a nonlinear kinematic constraint in spherical coordinates allowing a uniform radial displacement of the nodes at the void surface, u r ¼ cst and arbitrary displacement in the orthoradial directions (i.e. u f 1 and u f 2 are arbitrary where r, F 1 and F 2 denote the spherical coordinates).
The cases of a spherically "constrained" and a totally "unconstrained" void geometry are compared in the following for a large number of cycles. In the results that are shown subsequently, we stop the calculations at 50 cycles since after this point we observe strong localization of the deformation and the numerical solution diverges significantly. A more detailed discussion about this point is carried out in the Appendix B.
In Fig. 11, we show (a) the normalized average equivalent von Mises stress per cycle s eq /s 0 and (b) the porosity f per cycle as a function of the number of cycles N r for triaxiality X S ¼ 3, amplitude u/L ¼ 1% and Lode angle q ¼ 0 for the unconstrained void shape and the spherically constrained void shape cases. As we can observe in Fig. 11a, s eq /s 0 exhibits a maximum value for the unconstrained void shape (filled square on the curve indicated the position of this Fig. 9. Evolution of porosity f as a function of the average first principal strain at 5 cycles in the cases of (a) u/L ¼ 5%, C2, C3, C4, C5 represents the first, second, third, fourth and fifth cycle, respectively. Fig. 10. Evolution of the aspect ratio w 1 as a function of the average strain ε 1 at 5 cycles in the cases of (a) X S ¼ 2/3 (u/L ¼ 5%) and (b) X S ¼ 3 (u/L ¼ 1%) with q ¼ 0 . The notation C1, C2, C3, C4, C5 represents the first, second, third, fourth and fifth cycle, respectively. Due to the axisymmetric loading (q ¼ 0 ) w 2 ¼ 1 during the entire deformation process. maximum) whereas the spherically constrained void calculation exhibits a continuous hardening as a function of N r . In order to explain this rather interesting difference between the unconstrained and the constrained void shape calculations, we look into the porosity evolution in Fig. 11b. In this figure, we observe that the increase of porosity, albeit significant, cannot (by itself) explain the softening response observed in the unconstraint void shape calculation in Fig. 11a, since both the unconstrained void case as well as the spherically constrained void case predict almost the same evolution of the porosity as a function of the number of cycles N r , but only the first exhibits a maximum in the s eq /s 0 . (The point where the maximum in s eq /s 0 curve is observed is denoted with a filled square in Fig. 11b.) The fact that the spherically constrained case exhibits also porosity ratcheting can be partially explained by the presence of elasticity and the loading-unloading response of the unit-cell, as discussed in detail in Devaux et al. (1997). In fact, in the spherically constrained void shape case porosity ratcheting is even more pronounced than in the unconstrained void shape case. This obviously indicates that void shape effects are critical for the understanding of the cyclic response of such unit-cells. In this regard, in Fig. 12, we show the evolution (a) of the void aspect ratios w 1 and w 2 and (b) of the ellipsoidicity ratio ε l . While w 2 x 1 due to the axisymmetric loading along direction 1 (see Table 1 for q ¼ 0 ), w 1 evolves significantly as a function of the number of cycles N r . The constrained void shape curve is also shown for clarity (i.e., w 1 ¼ w 2 ¼ 1 during the deformation process). Perhaps more importantly, the ellipsoidicity ratio ε l dwhich serves to measure the deviation of the void shape from an ideal ellipsoidal shapedincreases significantly as a function of N r , in Fig. 12b. This indicates that not only the void becomes non-spherical but also deviates significantly from an ellipsoidal shape. Detailed contours of the underlying void geometry will be shown in the following. Thus, void geometry appears to be of crucial importance in the cyclic response of such unit-cells. This further suggests that small strains calculations may be insufficient for the analysis of cyclic loadings, at least in the present context where we attempt to analyze the effect of the underlying microstructure upon the average cyclic response of the unit-cell.

Effect of the loading and the initial void shape
In this section, we carry out a parametric study in order to investigate the influence of the loading conditions, i.e. the stress triaxiality and the Lode angle, as well as the effect of the initial void shape on the cyclic behavior of the periodic porous material. In the following results, the hardening exponent is set to N ¼ 10, the Young's modulus E ¼ 1000s 0 and the Poisson ratio n ¼ 0.3.   Fig. 12. Evolution of (a) the aspect ratios w 1 and w 2 and (b) the ellipsoidicity ratio ε l as a function of the number of cycles Nr for an unconstrained void and a spherically constrained void in the case of X S ¼ 3 and q ¼ 0 (u/L ¼ 1%). The filled square (-) on the graphs indicates the points where maximum equivalent stress is observed. less important (at least at the range of cycles considered here) but still non-negligible. In particular, the average stress s eq /s 0 exhibits no softening for the lower stress triaxiality X S ¼ 2/3 (Fig. 13a) contrary to the high stress triaxiality X S ¼3 (Fig. 13b). In addition, in Fig. 13b, for X S ¼ 3 and an unconstrained void shape (continuous line), we find that the s eq /s 0 response for q ¼ 0 depicts a more pronounced decrease than the two other cases, i.e., q ¼ 30 and q ¼ 60 . Again the spherically constrained void case shows no softening for any of the Lode angles considered here. These results indicate that primarily the stress triaxiality and secondary the Lode angle affect critically the average stress-strain response of the unitcell.

Effect of the stress triaxiality and the Lode angle
In order to address the role of the stress triaxiality and Lode angle upon the evaluation of s eq /s 0 , we show in Fig. 14 the porosity f as a function of the number of cycles N r for the same set of stress triaxialities and Lode angles for both the unconstrained void and the spherically constrained void shape. In Fig. 14a, for X S ¼ 2/3, we observe that even though the porosity evolves weakly (quantitatively) as a function of N r for all Lode angles considered, it tends to increase for the unconstrained void while decrease for the spherically constrained void, exhibiting a markedly different qualitative response in these two cases. On the other hand, for X S ¼ 3 in Fig. 14b we observe a very important increase of f (almost three times more that its initial value) but a less pronounced dependence on the Lode angle q. Moreover, due to the fact that the porosity evolution for X S ¼ 3 is quite similar for both unconstrained and constrained void shapes, we deduce again that the increase of porosity is not the only reason for the softening response observed in the unconstrained void shape calculation in Fig. 13b.
Next we examine the void geometry changes as a function of the stress triaxiality X S and the Lode angle q. More specifically, Fig. 15 shows the evolution of the aspect ratios w 1 and w 2 , as a function of the number of cycles N r for the same set of stress triaxialities and Lode angles considered previously. Obviously, there is no evolution of these quantities for spherically constrained void. The main observation in the context of Fig. 15 is that the evolution of the aspect ratios w 1 and w 2 becomes significant with increasing number of cycles. It is, in fact, observed that due to the applied finite deformations, the shape of the void changes from the very first cycle and it tends to grow further as the number of cycles increases. Even more interestingly, the largest change in the void shape occurs for higher stress triaxialities, i.e., X S ¼ 3, as shown in Fig. 15b, contrary to the case of X S ¼ 2/3 where both aspect ratios increase but in weaker manner. This result is not intuitive if one extrapolates the knowledge obtained in the context of monotonic loadings (see for instance Danas and Ponte Castañeda (2009a), Danas and Aravas (2012) and Srivastava and Needleman (2012)), where the largest void shape changes occur for lower stress triaxialities. At this point, we note that the evolution of the void shape does not describe adequately the deformation mechanisms near the void surface. In fact, for most of the computations presented here (except for the case of q ¼ 60 ) significant localization of the deformation occurs at the surface of the void. To illustrate this, we show, in Fig. 16 contours of the deformed unit-cell at 40 cycles for Lode angles q ¼ 0,30,60 and stress triaxiality (aec) X S ¼ 2/3 and (def) X S ¼ 3. To emphasize further the relative magnitude of each of the stress components js i j (i ¼ 1,2,3), we explicitly show them at the top corner of the displayed unit-cell. In particular, we observe a Fig. 13. Normalized maximum average equivalent von Mises stress evolution for unconstrained void/spherically constrained void in the case of (a) u/L ¼ 5%, X S ¼ 2/3 and (b) u/ L ¼ 1%, X S ¼ 3 as a function of the number of cycles Nr. The filled square (-) on the graphs indicates the point where maximum stress is observed. Fig. 14. Porosity evolution for unconstrained void/spherically constrained void in the case of (a) u/L ¼ 5%, X S ¼ 2/3 and (b) u/L ¼ 1%, X S ¼ 3 as a function of the number of cycles Nr.
The filled square (-) on the graphs indicates the point where maximum stress s eq /s 0 is observed. strong localization of the deformation (strains exceeding 60%) in a small zone of the void surface whose size depends upon the mesh size. While for q ¼ 0 , the localization strains lie in the plane 2À3 since the applied stress is axisymmetric along the 1-direction, for q ¼ 30 the localization zone is smaller (see for instance Fig. 16e) but still lying mainly in the plane 2À3. We should mention at this point that even though we have observed no maximum equivalent stress for the lower stress triaxiality X S ¼ 2/3 (see 13a), the deformation localization is more than present (Fig. 16a,b) and a critical event is expected to occur in these cases. Unfortunately, due to this strong localization of deformation the numerical calculation had to be stopped due to convergence issues as discussed further in the Appendix B, while a more appropriate formulation perhaps using non-local constitutive models (e.g., strain gradient plasticity models Forest et al. (2011), , Niordson (2013, 2014) or even more general non-local criteria Feld-Payet et al. (2011)) are needed in this case (see for instance Feld-Payet et al. (2011)). Those models should then to be combined with appropriate remeshing techniques but such a study is beyond the scope of the present work.
This same type of localization has been observed in all computations, i.e. for X S ¼ 2/3 and X S ¼ 3 except for q ¼ 60 . In this last case, Figs. 16c and f show that the void elongates significantly along the x 3 axis, i.e., along the direction of the minimum absolute stress component (since for q ¼ 60 , js 1 j ¼ js 2 j > js 3 j). It should be noted here that the observed localization affects only a small region of the void surface and inevitably leads to strong mesh dependence at the local level after localization occurs. As is detailed in the Appendix B, however, this mesh dependence affects local quantities (such as the aspect ratio and ellipsoidicity) and not average quantities such as the average stress and strain in the unit-cell, as well as the porosity evolution. It is worth noting at this point that recent experiments by continuous X-ray tomography (Hosokawa et al., 2012(Hosokawa et al., , 2013 have revealed the effect of stress triaxiality upon void shape and growth, albeit in monotonic loading conditions. As suggested by the present numerical calculations, the void shape changes under cyclic loadings are very different when compared to those obtained for monotonic loadings (Danas and Ponte Castañeda, 2009b;Srivastava and Needleman, 2012). Hence, a study similar to that of Hosokawa et al. (2012), but in cyclic loading conditions, could in fact shed light to the observed void shape effects.
To assess further the effects of the localization upon the void shape changes, we show in Fig. 17, the ellipsoidicity ratio ε l as a function of the number of cycles N r . We recall here that as ε l deviates from unity the void tends to diverge from an ideal ellipsoidal shape. In these graphs, the ellipsoidicity ratio reaches high values (ε l a1:5) for the cases X S ¼ 2/3, q ¼ 0 and X S ¼ 3, q ¼ 0, respectively, as a result of the corresponding deformation localization around the pore surface. The same remark can be made for X S ¼ 2/3 and q ¼ 30 but with somewhat lower values of the ellipsoidicity, in the order of ε l (1:3. In contrast, the ellipsoidicity is the smallest (and less than 1.2) for the combination q ¼ 30 and X S ¼ 2/3 as well as for q ¼ 30,60 and X S ¼ 3. This is in agreement with the contours of Fig. 16c,e,f, where the void shape does not deviate significantly from an ellipsoidal shape.

Effect of initial void shape
In this section, we investigate the effect of the initial void shape on the material cyclic response for initial porosity f 0 ¼ 0.01, triaxiality X S ¼ 3, amplitude u/L ¼ 1%, Lode angle q ¼ 0 ,60 and N ¼ 10.
In the following, only calculations with unconstrained void shapes will be shown. Fig. 18 shows evolution curves of (a) the average equivalent stress s eq /s 0 and (b) the porosity f as a function of the number of cycles for three different void geometries, i.e.
The initial void shapes that are different than a sphere have been chosen such that they remain axisymmetric during the deformation process, i.e., w 0 1 ¼ 1=3; w 0 2 ¼ 1 is a prolate void and w 0 1 ¼ 3; w 0 2 ¼ 1 is an oblate void, whose symmetry axis remains aligned with the maximum principal stress s 1 in this case. As observed in Fig. 18a, the average s eq /s 0 of the unit-cell exhibits a maximum only for the initially spherical void (w 0 1 ¼ w 0 2 ¼ 1) but not for the two other cases. Moreover, the unit-cell with w 0 1 ¼ 1=3; w 0 2 ¼ 1 is the stiffest of all three cases considered here. In Fig. 18b, the effect of the initial void shape upon the porosity evolution is also significant. It is worth emphasizing that the initially spherical void case (w 0 1 ¼ w 0 2 ¼ 1) shows the weaker porosity growth as a function of N r even though it is the only case that exhibits a maximum in the s eq /s 0 curve.
Similarly, Fig. 19 shows evolution curves of (a) the average equivalent stress s eq /s 0 and (b) the porosity f as a function of the number of cycles for three different values of the initial shape, i.e., w 0 1 ¼ w 0 2 ¼ 1; w 0 1 ¼ w 0 2 ¼ 1=3; w 0 1 ¼ w 0 2 ¼ 3 and for X S ¼ 3 and q ¼ 60 . The initial void shapes have also been chosen to remain axisymmetric during the entire deformation process, such as their symmetry axis remains aligned with the minimum principal stress s 3 in this case. For this configuration, in Fig. 19a, a maximum s eq /s 0 is observed for two out of three cases, i.e. for the initially spherical case (w 0 1 ¼ w 0 2 ¼ 1) and the initially prolate case (w 0 1 ¼ w 0 2 ¼ 3). In Fig. 15. Evolution of the aspect ratios w 1 (À), w 2 (À) for unconstrained void/spherically constrained void in the case of (a) u/L ¼ 5%, X S ¼ 2/3 and (b) u/L ¼ 1%, X S ¼ 3 as a function of the number of cycles Nr. The filled square (-) on the graphs indicates the point where maximum stress s eq /s 0 is observed. Fig. 19b the evolution of porosity f is less sensitive upon the initial void shape with that for w 0 1 ¼ w 0 2 ¼ 1 being the strongest among the three cases considered here contrary to the previous case of q ¼ 0 in Fig. 18.
In order to assess in a more visual and comprehensive way the evolution of the void shape geometry, Fig. 20 shows initial and final (at 50 cycles) void shapes for the same set of initial void geometries used in the previous figures of this section. The main observation in the context of this figure is that for stress triaxiality X S ¼ 3, the void tends to elongate in the direction of the minimum applied principal stress, i.e. js 2 j ¼ js 3 j for q ¼ 0 and js 3 j for q ¼ 60 , except in the case of q ¼ 60 and w 0 1 ¼ w 0 2 ¼ 1=3 (initially prolate case) where the void elongates along the direction of the maximum applied stress js 1 j ¼ js 2 j. This effect should not be confused with the similar effect observed for monotonic loadings and high stress triaxialities, since it is attributed to the fact that, in each cycle, the void shape evolves faster during the compressive step (X S 0) than the tensile step (X S ! 0), as already discussed in the context of Fig. 10 and is present also in lower stress triaxialities. At this point, it is worth noting that the direction of the void elongation is consistent (except for q ¼ 60 and w 0 1 ¼ w 0 2 ¼ 1=3) with a mode I loading of a crack. In other words, the voids tend to give a crack-type shape, as observed more clearly, in Fig. 20a,c and e, where a rather sharp crack tip starts to form in a direction perpendicular to the maximum applied average Fig. 16. Contours of the maximum principal logarithmic strain at 40 cycles in the case of X S ¼ 2/3 (u/L ¼ 5%) for (a) q ¼ 0 and (b) q ¼ 30 and (c) q ¼ 60 and in the case of X S ¼ 3 (u/ L ¼ 1%) for (d) q ¼ 0 and (e) q ¼ 30 and (f) q ¼ 60 . Fig. 17. Ellipsoidicity ratio ε l for an unconstrained void and a spherically constrained void in the case of (a) X S ¼ 2/3 (u/L ¼ 5%) and (b) X S ¼ 3 (u/L ¼ 1%) as a function of the number of cycles Nr. The filled square (-) on the graphs indicates the point where maximum equivalent stress is observed. stress in the unit-cell, i.e., as if they voids were under a mode I loading state. This observation simply indicates that the introduction of a pressure dependent defect (such as a pore) even if this defect is initially smooth (spherical or ellipsoidal), tends to diverge rapidly to a crack-type geometry and is expected to lead to the more common fatigue crack growth at larger number of cycles. Unfortunately at this level, severe mesh distortion did not allow us to continue further our calculation due to the strong localization of deformation. Remeshing of the final geometry together with nonlocal constitutive laws (Feld-Payet et al., 2011) might be required to proceed to higher number of cycles, however, such a calculation is beyond the scope of this work and is not pursued here.
At this point it is worth mentioning that while in monotonic loadings and under large triaxialities, the voids tend to grow rather uniformly and finally coalesce (see for instance Thomason (1985); Pardoen and Hutchinson (2000); Gologanu et al. (2001); Benzerga (2002)), in the present case of cyclic loading conditions coalescence of neighboring pores is not really due to void growth but rather due to crack initiation and propagation, which is also present at lower stress triaxialities. As observed in Figs. 16 and 20, the voids do not grow uniformly but rather a crack is created in a very thin localized zone for both stress triaxialities analyzed in this work. This crack is then expected to propagate and coalesce with the neighboring pores but such an analysis is highly mesh dependent and not pursued further here as discussed previously. As a result in the context of cyclic loadings, the localization of plastic strain does not take place within horizontal ligaments spanning the size of the pore and the spacing of the neighboring pores as in the usual studies of coalescence in monotonic loadings (Pardoen and Hutchinson, 2000) but rather in a confined thin region around a small portion of the void surface.

Preliminary results on combined isotropic-kinematic hardening
At this point, it would be constructive to for a matrix material following combined isotropic and nonlinear kinematic hardening. The range of parameters needed to investigate in a full extent the effect of kinematic hardening upon the cyclic response of the unitcell is large and thus prohibitive for the present work and will be presented elsewhere in a subsequent study. Nonetheless, as we will see below the introduction of nonlinear kinematic hardening in the matrix phase does not alter the qualitative character of the deformation zone and localization around the void but only the quantitative aspects.
More specifically, let us consider the same isotropic hardening presented in equation (17) together with a nonlinear kinematic hardening with two-backstresses X (1) and X (2) (c.f. Nouailhas et al., 1985), i.e., In this expression, D p M is the second-order plastic strain-rate tensor in the matrix phase (the subscript M has been used to explicitly distinguish between strain-rate in the matrix phase and in the case of X S ¼ 3 (u/L ¼ 1%) and q ¼ 0 as a function of the number of cycles Nr. The filled square (-) on the graphs indicates the point where maximum equivalent stress is observed. Fig. 19. (a) Normalized maximum average equivalent von Mises stress evolution for several initial void shapes (w 0 Porosity evolution for several initial void shape (w 0 the average strain-rate in the unit-cell) and ε p is the accumulated plastic strain defined in the context of equation (17). Furthermore, motivated by the choices made in Besson and Guillemer-Neel (2003), we set the values C 1 ¼ 2.5 Â 10 À2 s 0 , C 2 ¼ 10C 1 and g 2 ¼ 10. These values correspond to a rather week kinematic hardening at the level of the matrix phase but as will be shown next they already lead to a strong quantitative (but not qualitative) effect upon the cyclic response of the unit-cell. The q h serves to parametrize the relative effect of kinematic hardening with respect to the isotropic hardening and takes here three values q h ¼ 0,0.5,1 (with q h ¼ 0 corresponding to the case of no kinematic hardening). As a representative result in the context of combined isotropic/ kinematic hardening, we have chosen the case of high stress triaxiality X S ¼ 3 (u/L ¼ 1%) and Lode angle q ¼ 0 . In Fig. 21, we show contours of the maximum principal strain for (a) q h ¼ 0 at N r ¼ 40 cycles (no kinematic hardening), (b) q h ¼ 0.5 at N r ¼ 60 cycles and (c) q h ¼ 1 at N r ¼ 100 cycles. The main observation in the context of this figure is that the deformation localization locus is the same with or without kinematic hardening, i.e., the strain localizes in the same position at the void surface in an attempt to create a mode-I crack-type geometry, as discussed in the previous sections. The difference however is quantitative, i.e., with the addition of kinematic hardening this deformation localization occurs at much larger number of cycles as shown in Fig. 21. In other words, the addition of a rather small value of kinematic hardening does not alter the character of the deformation localization in the zone around the void, but rather changes the number of cycles needed to reach this localized state. In a sense, it "decelerates" the localization procedure, thus having a beneficiary effect upon the cyclic response of the material.
It should be further pointed out that similar effects have been found for other values of the stress triaxiality and the Lode parameter but are not shown here for brevity. These preliminary results call for an in-depth quantitative analysis of the effect of kinematic hardening upon the cyclic response of such unit-cells and will be presented in a future work.

Discussion
In this work we have investigated the effects of cyclic loading conditions upon microstructure evolution and material softening/ hardening using a cubic periodic unit-cell comprising a single spherical (or ellipsoidal) void at the center with volume fraction 1% that is subjected to triaxial finite deformations such that the stress triaxiality and the Lode angle are kept constant during each step of the cycle. It has been found that the void shape changes are asymmetric when subjected to positive and negative stress triaxialities. This is exhibited by a permanent irreversible void shape change from the very first cycle, which in turn leads to porosity ratcheting and an asymmetric average stress-strain response. It Fig. 20. Cross-sections of undeformed (continuous line) and deformed (dashed lines) void geometries at 50 cycles and stress triaxiality X S ¼ 3 in the case of q ¼ 0 with initial void shape (a) w 0 1 ¼ w 0 2 ¼ 1 and (b) w 0 1 ¼ 1=3; w 0 2 ¼ 1 and (c) w 0 1 ¼ 3; w 0 2 ¼ 1 in the case of q ¼ 60 with initial void shape (d) w 0 should be noted that this asymmetry does not lead to a Bauschinger effect but is strongly related to the evolution of the underlying void volume and shape. Furthermore, we have observed that both the stress triaxiality and the Lode angle can have strong effects upon the cyclic response of the unit cell. It has been found that for initially spherical voids the average stress in the unit-cell exhibits a maximum as a result of the critical void shape changes at high stress triaxialities, but not due to porosity ratcheting alone. To establish this conclusion, we have carried out additional calculations on the same unit-cell but constraining the shape of the void to remain spherical during the entire deformation process. In this case, we have also obtained significant porosity ratcheting (very similar and even larger than that of the unconstraint void), but no maximum was observed in the average stress response of the unitcell. This suggests that void shape changes are critical in the cyclic response of the unit-cell. In addition, the observed void shape changes have led to significant localization of the deformation near the void surface. Strains larger than 50% have been observed near the void surface even when the average strain amplitude in the unit-cell has been of the order of 1%. This result further revealed the importance of carrying out the cyclic analysis under a finite deformation framework. Moreover, we have investigated the effect of initially nonspherical (i.e., ellipsoidal) voids on the cyclic response of the unit-cell. The main outcome of this parametric analysis has been that the void tends to elongate in the direction of the minimum absolute applied stress. This observation has been different only in a few distinct cases where the void shape was positioned in a significantly non-optimal direction with respect to a mode I loading direction. This void elongation, in turn, leads to a crack-shape microstructure subjected to mode I cyclic loading state and to eventual coalescence for both lower and higher stress triaxialities. Nevertheless, the strain fields in that late stage of coalescence are very different from those observed in the context of monotonic loading states and high triaxialities. In the context of cyclic loads, the localization zone is confined in a very thin region around a small portion of the void surface.
Finally, as discussed before, the cyclic behavior of materials exhibits different stages where isotropic and/or kinematic hardening (Nouailhas et al., 1985;Chaboche, 2008) are the main cyclic mechanisms. The present study provides a partial perspective upon the cyclic response of a periodic porous material by considering mainly isotropic hardening and only a test case with combined nonlinear isotropic and kinematic hardening. A complete study that analyzes the relative importance of isotropic and kinematic hardening upon the cyclic response of the unit-cell should be considered further. One should also mention that anisotropy of the matrix phase (e.g., crystal plasticity) is also expected to have an effect upon the cyclic response of such porous unit-cells. Nonetheless, the above-presented preliminary calculations show that, for instance, the presence of kinematic hardening exhibits similar qualitative behavior with the case of no kinematic hardening but tends to decelerate the initiation of localization of deformation around the void surface. Such work is in progress and will be reported elsewhere. X S ¼ 3, amplitude u/L ¼ 1%, Lode angle q ¼ 60 , and a set of hardening exponents N ¼ 5,10,∞ (where the last one corresponds to the perfect plasticity case). Initially, we discuss the effect of the hardening exponent for a realistic value of the elastic moduli, e.g., Young's modulus, E~1000s 0 (for instance s 0~2 00 MPa corresponding to steel) and then we also present a brief result of porosity ratcheting for E/s 0 ¼ 50,000 in order to address the observations made by Devaux et al. (1997) in the limit of vanishing elasticity (i.e., E/s 0 / ∞).
As we can observe in Fig. 22a, the average equivalent von Mises stress per cycle s eq /s 0 decreases monotonically with respect to N. In particular, for the case of the unconstrained void shape, it exhibits a maximum value for s eq /s 0 for N ¼ 10 but not for N ¼ 5, at least up to 50 cycles where the calculations are terminated. On the other hand, as expected, s eq /s 0 shows no maximum for the N ¼ 5 and N ¼ 10 for the spherically constrained shape. For the case of perfect plasticity, Fig. 21. Contours of the maximum principal logarithmic strain in the case of X S ¼ 3 (u/L ¼ 1%) and q ¼ 0 for (a) q h ¼ 0 at N r ¼ 40 cycles (no kinematic hardening), (b) q h ¼ 0.5 at N r ¼ 60 cycles and (c) q h ¼ 1 at N r ¼ 100 cycles.
i.e. N / ∞, both the unconstrained and the constrained void shape exhibit a softening from the very first cycle. In Fig. 22b, however, the evolution of the porosity exhibits nonmonotonic dependence upon N as a function of the number of cycles N r . The evolution of f is higher for N ¼ 10 and lower for N ¼ 5, with the one corresponding to N ¼ ∞ intersecting the other two curves at about 10 cycles. As we will see below this intersection could be attributed to the localization of deformation which is more pronounced in the case of N ¼ ∞, but this point should further be studied in the future.
In Fig. 23, we show evolution curves for (a) the aspect ratios w 1 ¼ w 2 (for q ¼ 60 ) and (b) for the ellipsoidicity as a function of the number of cycles N r for the same set of hardening exponents used before (N ¼ 5,10,∞). The evolution of the aspect ratios in Fig. 23a are shown to be almost independent of N as N ! 10, while they tend to evolve slower for larger N ¼ 5. In contrast the evolution of ellipsoidicity in Fig. 23b exhibits a non-monotonic dependence on N especially at large number of cycles taking the largest values for N ¼ 10.
Finally, in Fig. 24, we consider the case of vanishing elasticity, i.e., E/s 0 ¼ 50,000 and record porosity ratcheting under the same loading state and the same set of hardening exponents (N ¼ 5, 10,∞). One observes that all results confirm (up to a given numerical accuracy since the Young's modulus is not infinity) the results by Devaux et al. (1997), i.e., there is no porosity ratcheting effect in the absence of elasticity. Rather interestingly, porosity ratcheting also becomes negligible for N ¼ 5,10. But this is related to the rather small macroscopic strain amplitude considered here as well as to the specific cyclic loading conditions, i.e., cyclic loads around zero average straining. Thus, if one compares these results with the previous ones (which correspond to E ¼ 1000s 0 ) then the effect of elasticity is critical for porosity ratcheting.

Appendix B. Mesh dependence
In this section, we discuss the dependence of mesh size upon the cyclic response of the periodic unit-cell. As we have observed on several calculations with q ¼ 0,30 (see Fig. 16), a localization of deformation appears in a small region of the void surface and thus, leads to strong mesh dependence at this zone. This led to nonconvergence of the elasto-plastic computations beyond 50 cycles. In Fig. 25, we show representative results for a loading u/L ¼ 1%, X S ¼ 3, q ¼ 30 and for three different mesh sizes at the pore surface, i.e., 0.052a, 0.037a and 0.026a (where a is the radius of the initially spherical pore), corresponding to 5 Â 10 4 , 15 Â 10 4 and 32 Â 10 4 degrees of freedom (dof), respectively. Fig. 25a shows the evolution of the maximum average equivalent von Mises stress as a function of the number of cycles N r for the same loading and the aforementioned three different mesh sizes. It is evident from Fig. 25a,b that the evolution of the average stress s eq /s 0 (as well as the position where maximum is attained) and the porosity f (which is a direct measure of the hydrostatic average plastic strain) as a function of N r is not sensitive to the mesh size. On the other hand, in Fig. 25c,d we observe that the aspect ratios w 1 and w 2 are much more sensitive to the mesh size. In particular, even though the evolution of w 1 and w 2 tend to converge for N r < 20, the curves diverge after this point as a result of the strong deformation localization observed in those cases. In addition, it is interesting to note that there is a monotonic decrease of the aspect ratios with increasing mesh size. Consequently, we can conclude that the mesh dependence affects more the local quantities (such as the aspect ratio) but not average quantities such as the average stress or strain, and the porosity evolution, which depends on the hydrostatic part of the average strain in the unit-cell.   23. (a) Evolution of the aspect ratios w 1 (À), w 2 (À) for several values of hardening exponent (N ¼ 5,10,∞) in the case of u/L ¼ 1%, X S ¼ 3, q ¼ 60 and (b) Ellipsoidicity ratio for several values of hardening exponent (N ¼ 5,10,∞) in the case of u/L ¼ 1%, X S ¼ 3, q ¼ 60 as a function of the number of cycles Nr.  s eq /s 0 , (b) porosity f, (c) aspect ratio w 1 and (d) aspect ratio w 2 for 5 Â 10 4 , 10 5 and 3 Â 10 5 degrees of freedom (dof) in the case of X S ¼ 3 and q ¼ 30 as a function of the number of cycles Nr.