A model for creep of porous crystals with cubic symmetry

Accepted Manuscript A model for creep of porous crystals with cubic symmetry A. Srivastava, B. Revil-Baudard, O. Cazacu, A. Needleman PII: DOI: Reference: S0020-7683(17)30056-2 10.1016/j.ijsolstr.2017.02.002 SAS 9460 To appear in: International Journal of Solids and Structures Received date: Revised date: Accepted date: 25 July 2016 16 December 2016 3 February 2017 Please cite this article as: A. Srivastava, B. Revil-Baudard, O. Cazacu, A. Needleman, A model for creep of porous crystals with cubic symmetry, International Journal of Solids and Structures (2017), doi: 10.1016/j.ijsolstr.2017.02.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. ACCEPTED MANUSCRIPT A model for creep of porous crystals with cubic symmetry A. Srivastavaa , B. Revil-Baudardb , O. Cazacub , A. Needlemana a Department of Materials Science & Engineering, Texas A & M University, College Station, TX, USA Department of Mechanical and Aerospace Engineering, University of Florida/REEF, Shalimar, FL, USA CR IP T b Abstract PT ED M AN US A model for description of the creep response of porous cubic single crystal is presented. The plastic potential is obtained by specializing the orthotropic potential of Stewart and Cazacu (Int. J. Solids Struct., 48, 357, 2011) to cubic symmetry. The crystal matrix material response is characterized by power law creep. The predictions of this porous plastic constitutive relation are presented for various values of stress triaxiality (mean normal stress divided by Mises effective stress) and various values of the Lode parameter L (a measure of the influence of the third invariant of the stress deviator). A strong influence of crystal orientation on the evolution of the creep strain and the porosity is predicted. For loadings along the < 100 > directions of the cubic crystal, void growth is not influenced by the value of the Lode parameter. However, for loadings such that the maximum principal stress is aligned with the [110] direction there is a strong influence of the values of the Lode parameter and the fastest rate of void growth occurs for shear loadings (one of the principal values of the applied stress deviator is zero). For loadings such that the maximum applied stress is along the [111] crystal direction the fastest rate of void growth corresponds to L=-1, while the slowest rate corresponds to L=1. These predictions are compared with corresponding predictions of the three dimensional finite deformation unit cell analysis of Srivastava and Needleman (Mech. Mater., 90, 10, 2015). It is found that the phenomenological model predicts the same trends as the cell model calculations and, in some cases, gives good quantitative agreement. CE Keywords: Porous material; Creep; Single Crystal; Anisotropy; Lode parameter 1. Introduction AC As is well-appreciated, porosity nucleation, growth and coalescence is the main mechanism of ductile fracture in structural metals. Porosity evolution can also play a role in determining the deformation response of structural metals and alloys in circumstances where fracture is not a main concern but the overall deformation is, for example in deformation processing of sintered materials. Therefore, a basic understanding of the evolution of porosity and its effect on the overall mechanical response is of widespread interest. In a wide range of circumstances, the voids of interest are of a size (say several microns and larger) where the surrounding material can be appropriately characterized by a continuum plasticity constitutive description. Preprint submitted to Elsevier February 4, 2017 ACCEPTED MANUSCRIPT AC CE PT ED M AN US CR IP T There is a long history of continuum based unit cell model calculations of void containing solids as well as of the development of phenomenological theories of porous plastic solids, see Tvergaard (1990), Needleman et al. (1992) and Benzerga and Leblond (2010) for reviews with a focus on ductile fracture applications. Most of the constitutive relations that have been developed for porous plastic solids have presumed that the surrounding matrix material can be regarded as isotropic. There are circumstances, however, where anisotropy of the matrix material may play a significant role, for example, for strongly textured polycrystalline solids and where the void is embedded in a single crystal or is surrounded by a few crystals. Using rigorous limit-analysis theorems, Benzerga and Besson (2001) derived an analytic yield function for a porous material containing randomly distributed spherical voids in a matrix obeying Hill (1948) orthotropic criterion. Within the same framework, Monchiet et al. (2008) studied the case of ellipsoidal voids and derived a closed-form orthotropic plastic potential. The combined effects of anisotropy and tension-compression asymmetry induced by twinning or non-Schmid effects on the dilatational response of porous textured polycrystals was investigated by Stewart and Cazacu (2011), Cazacu and Stewart (2013). These authors analytically solved a limit-analysis problem for the cases of spherical and cylindrical void geometries, respectively, and developed appropriate orthotropic plastic potentials. However, for the case of porous single crystals analytical derivation of a plastic potential poses challenges that, at least to now, have been insurmountable. If the plastic deformation of the fully-dense single crystal is described using either a Bishop and Hill (1951) type model or the regularized form proposed by Arminjon (1991) with the exponent n 6= 2, the plastic dissipation cannot be expressed in closed-form. Therefore, it is impossible to solve the limit-analysis problem analytically and derive a closed-form expression for the plastic potential. This was stated in Paux et al. (2015), who proposed an ad-hoc modification of the Gurson (1975) isotropic yield function. An alternative approach to model the mechanical response of porous materials is based on the homogenization method developed for non-linear composites by Ponte-Castaneda (see for example, Ponte Castaneda (2002)). This method is based on the equivalence response of the solid under consideration with a linear-comparison composite solid described by a potential quadratic in stresses. It was applied by Idiart and Ponte Castaneda (2007) for the study of porous single crystals containing cylindrical voids subject to anti-plane loadings, and more recently by Mbiakop et al. (2015) for two-dimensional plane strain loadings. As pointed out by Mbiakop et al. (2015), for anisotropic crystal plasticity based on a powertype law description with exponent n 6= 2, no analytic solution exists even for hydrostatic loadings. Nevertheless, Mbiakop et al. (2015) were successful in obtaining numerical plastic potential surfaces for various loadings. However, no results were reported for the time evolution of plastic strain or porosity under creep loading. For the case of a porous single-crystal with a matrix obeying a quadratic (i.e. n = 2) Bishop and Hill (1951) relation, Han et al. (2013) used the linear-comparison composite solid method to obtain an approximate analytical plastic potential that is quadratic in the components of stress. Han et al. (2013) also compared the model predictions with finite-element cell calculations for different crystal orientations. Three-dimensional cell model calculations exploring the effect of crystal induced anisotropy on the stress state dependence of porosity evolution were reported in Wan et al. 2 ACCEPTED MANUSCRIPT AC CE PT ED M AN US CR IP T (2005); Yu et al. (2010); Ha and Kim (2010); Yerra et al. (2010); Lebensohn and Cazacu (2012); Han et al. (2013); Srivastava and Needleman (2012, 2013, 2015). Finite element modeling of the plastic deformation of single crystals for example fcc crystals requires accounting for slip on each of the twelve available slip systems. This additional computational complexity limits the use of such a model in applications. In addition, although for a rate independent single crystal obeying a Schmid slip system relation the yield surface (and hence the flow potential) is faceted and has sharp corners, rate dependence rounds off the corners and gives rise to a smoother flow potential surface when multiple slip systems are significantly active, Rice (1970). The aim of this paper is to provide a simple phenomenological model for representing the creep response of porous cubic single-crystals. In particular, for simplicity and to keep the expressions close to those derived analytically from limit analysis, we account for crystal anisotropy but not for the discreteness of slip systems. The proposed phenomenological model is obtained by specializing the orthotropic potential derived by Stewart and Cazacu (2011) to the case of cubic symmetry. To account for rate-effects, we use the approach proposed by Pan et al. (1983). We compare the predictions of the proposed phenomenological model with the three dimensional single crystal unit cell results of Srivastava and Needleman (2015). Their finite deformation finite element calculations were carried out for an fcc single crystal containing a single initially spherical void. The deformation of the matrix was modeled by a crystal plasticity (Asaro and Needleman, 1985) framework with a power law viscous creep relation for the matrix material. The unit cell was subject to creep loading, i.e. a fixed stress state, for a range of values of the imposed stress triaxiality, the ratio of the first to second stress invariants, and a range of imposed values of the Lode parameter, a measure of the third stress invariant. The results of Srivastava and Needleman (2015) showed a strong effect of anisotropy and stress state on the evolution of the overall creep strain and porosity. As expected, the predicted response was found to be sensitive to the value of the applied stress triaxiality. For the [100] crystal orientation that gives rise to nearly isotropic response, no effect of the Lode parameter on the dilatational response was observed. On the other hand, for anisotropic crystal orientations, a significant influence of the Lode parameter was found on the creep response of the porous crystals even at a high value of the stress triaxiality. In this paper, using the proposed phenomenological model the effect of crystal orientation is analyzed for various creep loading conditions. The model predictions for the overall creep response and for porosity evolution are compared with the corresponding cell model results of Srivastava and Needleman (2015). Our analytical results show a strong effect of crystal orientation and imposed stress state on the evolution of overall creep strain and of porosity. The two important distinction between the cell model calculations and the simple phenomenological model are: (i) in the cell model calculations of Srivastava and Needleman (2015) the orientations of the slip systems evolve, whereas in the results based on the plastic potential of Stewart and Cazacu (2011) the anisotropy is fixed throughout the deformation history; (ii) the cell model calculations account for void-void interactions, whereas in the simple phenomenological model any such interactions are ignored. Nevertheless, key features of the phenomenological predictions are consistent with those obtained from the cell model 3 ACCEPTED MANUSCRIPT calculations. 2. Formulation CR IP T To describe the creep response of porous cubic single-crystals, we specialize the orthotropic plastic flow potential of Stewart and Cazacu (2011) to cubic symmetry. The plastic potential of Stewart and Cazacu (2011) is briefly described in Section 2.1 and a model for creep of porous crystals with cubic symmetry is proposed in Section 2.2. φ (σij , f ) = m̂ 2 AN US 2.1. The plastic potential for orthotropic porous solids of Stewart and Cazacu (2011) Stewart and Cazacu (2011) used a kinematic limit analysis approach in conjunction with the Hill-Mandel lemma (Hill, 1967; Mandel, 1972) to derive an analytical expression for the plastic potential of an orthotropic rate independent plastic solid containing randomly distributed spherical voids. The plastic behavior of the matrix (void-free solid) was taken to be governed by a relation that accounts for plastic tension-compression asymmetry but is pressure-insensitive, (Cazacu et al. (2006)). The Stewart and Cazacu (2011) plastic potential has the form 2 3  X |σ̃i | − kσ̃i i=1 σxT + 2f cosh  3σm hσxT  − (1 + f 2 ) (1) ED M where f is the void volume fraction (or porosity), k is a material parameter accounting for the tension-compression asymmetry in plastic deformation, σxT is the uniaxial tensile yield strength along an axis of orthotropy, σm = tr (σij ) /3 and σij are Cartesian components of the Cauchy stress tensor. In Eq. (1), σ˜1 , σ˜2 , σ˜3 are the principal values of the transformed stress tensor 0 σ̃ij = κijkl σkl (2) AC CE PT where κijkl is a fourth rank symmetric orthotropic tensor and σij0 is the deviator of the Cauchy stress tensor, σij , i.e. σij0 = σij − δij σm . In the coordinate system with axes along the orthotropy directions (e.g. for a rolled sheet, these axes are the rolling, transverse, and through thickness directions), the fourth rank tensor κ in Voigt notation is written as   κ11 κ12 κ13 0 0 0  κ12 κ22 κ23 0 0 0     κ13 κ23 κ33 0 0 0    κ= (3)  0 0 0 κ 0 0 44    0 0 0 0 κ55 0  0 0 0 0 0 κ66 Also, in Eq. (1), m̂ is a material parameter and is expressed in terms of the components of the anisotropy tensor, κ, and the tension-compression asymmetry parameter, k, as 1 m̂ = q (|Φ1 | − kΦ1 )2 + (|Φ2 | − kΦ2 )2 + (|Φ3 | − kΦ3 )2 4 (4) ACCEPTED MANUSCRIPT Here, Φ1 = (2κ11 − κ12 − κ13 ) /3, Φ2 = (2κ12 − κ22 − κ23 ) /3 and Φ3 = (2κ13 − κ23 − κ33 ) /3. The material parameter h in Eq. (1) depends on both the anisotropy coefficients and the sign of the mean stress σm and is given by r r h= (4t1 + 6t2 ) (5) 5 r=    1 m̂2 1 m̂2 3 3k2 −2k+3 3 3k2 +2k+3

CR IP T with if σm < 0 if σm ≥ 0 2 2 2 t1 = 3 B13 B23 + B12 B23 + B12 B13 + 2B12 + 2B13 + 2B23 2 2 2 t2 = B44 + B55 + B66
(6) (7) (8) AN US where Bij are the components of inverse of the deviator of the fourth rank tensor κ.1 Thus, plastic flow and hence porosity evolution depends on orientation and whether the applied loading is tensile or compressive. For the matrix material, i.e., for f = 0, Eq. (1) reduces to the plastic potential of Cazacu et al. (2006). CE PT ED M 2.2. Model for creep of porous cubic crystals To describe the plastic flow rule for a creeping cubic single-crystal in the presence of voids, we specialize the plastic flow potential of Stewart and Cazacu (2011) presented in Section 2.1. We presume that the tension-compression asymmetry is negligible, so that k = 0 in Eq. (1). The potential in Eq. (1) was derived presuming a non-hardening rate independent matrix material. We explore the extent to which the form of such a plastic flow potential can be used to describe highly rate dependent deformation. If it can, it means that the forms of plastic flow potentials suggested by such analyses can, at least in some circumstances, provide a basis for modeling material behavior well outside the range of matrix material responses for which they were derived. Denoting the base vectors of the Cartesian laboratory frame by Ci and the base vectors of the (Cartesian) crystal frame by ci where for the cubic crystals considered c1 = [001], c2 = [010] and c3 = [001], (9) Ci = Qij cj ci = QTij Cj AC where Q is a rotation matrix. The stress tensor Σ can be written as Σ = Σij Ci ⊗ Cj = σij ci ⊗ cj (10) 1 See Cazacu et al. (2006, 2010); Stewart and Cazacu (2011) for expressions relating the components of B to the components of the tensor κ. 5 ACCEPTED MANUSCRIPT On the crystal the stress components are specified on the laboratory axes Σij and we need to obtain the stress components on the crystal axes. From Eqs. (9) and (10) Σij Ci Cj = Σij Qim cm Qjn cn = σmn cm cn (11) σmn = Qim Qjn Σij (12) so that CR IP T The rotation matrix Q simply depends on the orientation of the crystal with respect to the laboratory/loading axes, Fig. 1. The flow potential for the porous crystal takes the form   σkk σ̃e2 − (1 + f 2 ) = 0 (13) φ = 2 + 2f cosh σM hσM AN US Note that σkk = Σkk and σM is an internal variable representing the average matrix flow strength, and σ̃e2 = m̂2 σ̃ij σ̃ij (14) M the deviatoric part of the transformed stress tensor, σ̃ij0 is related to the deviatoric part of the Cauchy stress tensor, σij , through Eq. (2). Also the Cauchy stress tensor, σij , in the crystal frame is related to the stresses specified in the laboratory frame, Σij , through Eq. (12). In the crystal frame, the anisotropic tensor κ satisfies κ11 = κ22 = κ33 , κ12 = κ13 = κ23 , κ44 = κ55 = κ66 ED so that (15) (16) (17) (18) σ̃12 = κ44 σ12 , σ̃13 = κ44 σ13 , σ̃23 = κ44 σ23 (19) CE PT 0 0 0 σ̃11 = κ11 σ11 + κ12 (σ22 + σ33 ) 0 0 0 σ̃22 = κ11 σ22 + κ12 (σ11 + σ33 ) 0 0 0 σ̃33 = κ11 σ33 + κ12 (σ11 + σ22 ) AC From Eqs. (16), (17) and (18), it follows that σ̃ij is traceless. Also, from the definitions in Section 2.1, Φ2 = Φ3 = −Φ1 /2 and so that and 2(κ11 − κ12 ) 3 (20) 2 3 1 3 β2 = = 3Φ21 2 (κ11 − κ12 )2 2 κ244 (21) Φ1 = m̂2 = v " s   2 # s    u 2 u 3 8 9Φ1 1 κ11 − κ12 1 1 t h= + 2 = 8 + 12 = 8 + 12 2 5 3 κ44 5 κ44 5 β 6 (22) ACCEPTED MANUSCRIPT with 2 β =  κ44 κ11 − κ12 2 (23) As shown in the Appendix, the expression for σ̃e2 , Eq. (14), reduces to σ̃e2 =

3  02 02 02 02 02 02 + σ23 + σ13 + 2β 2 σ12 + σ33 σ11 + σ22 2 (24) AN US CR IP T with β as defined in Eq. (23). The anisotropy coefficients κij characterize the anisotropy of the plastic response of the matrix material (the fully-dense crystal). Since σ̃e2 is homogeneous of degree zero in the coefficients κij , the response is the same if the coefficients κij are replaced by aκij , where a is any positive constant. Therefore, without loss of generality one of the coefficients κij can be set equal to unity; for example we can take κ11 = 1. Moreover, the plastic response depends only on the ratio of κ12 to κ44 , or alternatively on β, see Eq.(23). For isotropy κ12 = 0 and κ44 = 1 so that m̂2 = 3/2 and h = 2 recovering the Gurson (1975) potential. As used by Pan et al. (1983) to extend the Gurson (1975) model to the case of viscoplastic behavior, rate dependent matrix material response is accounted for here by incorporating a rate dependent flow strength σM in Eq. (13). Thus, for power law creep, the matrix strain rate is  n σM ˙M = η (25) σ0 PT ED M where (˙) denotes differentiation with respect to time and, η, σ0 and n are material parameters. Given the orientation matrix Q and the anisotropy coefficients κ11 , κ12 , κ44 , together with the porosity f and the stress components Σij all quantities in the flow potential Eq. (13) are explicitly determined. Neglecting elasticity, the components of the rate of deformation tensor D on the laboratory axes are obtained via Dij = Λ̇ ∂φ ∂Σij (26) AC CE where Λ̇ is a scalar parameter with dimension 1/time. The evolution of the void volume fraction is obtained from conservation of mass with an incompressible matrix so that f˙ = (1 − f )Dkk (27) To identify Λ̇ we use the equivalence of overall dissipation W and matrix plastic dissipation, i.e. ∂φ W = Σij Dij = (1 − f )σM ˙M = Σmn Λ̇ (28) ∂Σmn The plastic dissipation W is non-negative for f ≤ 1. Hence, W (1 − f )σM ˙M (1 − f )ησM (σM /σ0 )n Λ̇ = = = (29) Σmn ∂Σ∂φmn Σmn ∂Σ∂φmn Σmn ∂Σ∂φmn 7 ACCEPTED MANUSCRIPT where the last expression presumes the power law creep relation, Eq. (25), for the crystal matrix. To obtain explicit expressions for Dij and f˙, we use   σ̃e2 Σkk Σkk ∂φ = 2 2 + 2f sinh (30) Σmn ∂Σmn σM hσM hσM
∂φ 2m̂2 2f = 2 κklmn Qpm Qqn Σ0pq (κklrs Qir Qis ) + sinh ∂Σij σM hσM Then, from Eqs. (26) and (29), Dij is given by Dij =  (1 − f )ησM (σM /σ0 )n ∂φ   2 Σkk Σkk ∂Σij 2 σσ̃2e + 2f hσ sinh hσM M M Σkk hσM  CR IP T and δij (31) (32) AN US with ∂φ/∂Σij given by Eq. (31) and from Eq. (27) the evolution of porosity is governed by   (1 − f )ησM (σM /σ0 )n Σkk 2f ˙   f = (1 − f ) 2 sinh (33) Σkk Σkk hσM hσM 2 σσ̃2e + 2f hσ sinh hσM M M 3. Numerical Results ED M The value of σM is obtained from the consistency condition, φ = 0, in Eq. (13). Note that the crystal orientation matrix Q only enters Eqs. (32) and (33) through the expression κklmn Qpm Qqn . Thus, it is straightforward to calculate Dij and f˙ for any crystal orientation. PT Table 1: Values of the stress triaxiality χ, the Lode parameter L, and the initial overall stresses Σi . AC CE χ 3 3 3 2/3 2/3 2/3 L Σ1 (MPa) Σ2 (MPa) Σ3 (MPa) -1.00 2750.00 2000.00 2000.00 0.00 2683.01 2250.00 1816.99 1.00 2500.00 2500.00 1750.00 -1.00 1000.00 250.00 250.00 0.00 933.01 500.00 66.99 1.00 750.00 750.00 0.00 The predictions of the constitutive model presented in Section 2.2 are explored for crystal orientations and stress states for which finite element finite deformation unit cell analyses were carried out in Srivastava and Needleman (2015). These predictions are then compared with the cell model results in Srivastava and Needleman (2015). The five crystal orientations are shown in Fig. 1. The primary orientations have the [100], [110] and [111] directions parallel to the main loading direction (the x1 -axis). Two 8 M AN US CR IP T ACCEPTED MANUSCRIPT PT ED Figure 1: Relative orientations of the cubic crystal directions (crystal frame) and the coordinate axes (laboratory/loading frame) for the five cases analyzed. The main loading direction is taken to be parallel to the x1 axis. AC CE secondary orientations for the [110] and [111] orientations are considered as shown in Fig. 1. As in Srivastava and Needleman (2015), the stress states analyzed are defined in terms of the values of the imposed stress triaxiality, χ, and the Lode parameter, L, which are given by 2Σ2 − Σ1 − Σ3 Σh L= (34) χ= Σe Σ1 − Σ3 with 1 p 1 Σe = √ (Σ1 − Σ2 )2 + (Σ2 − Σ3 )2 + (Σ3 − Σ1 )2 Σh = (Σ1 + Σ2 + Σ3 ) (35) 3 2 In Eq. (34), for the calculation of the Lode parameter L the principal stresses are ordered such that Σ1 ≥ Σ2 ≥ Σ3 . The stress states for which creep calculations are carried out for the constitutive model in Section 2.2 are given in Table 1. The values of Σi are held fixed, as is their directions along the laboratory coordinate axes. The evolution of creep strain and porosity is calculated using 9 ACCEPTED MANUSCRIPT CR IP T Eqs. (32) to (33). Attention is restricted to triaxiality values χ = 3 and χ = 2/3. Srivastava and Needleman (2015) also carried out calculations for χ = 1/3 but it is well-know that potentials having the form in Eqs. (13) and (1) are not accurate for low values of χ where void shape effects play a significant role, hence χ = 1/3 is not considered here. The crystal orientation with respect to loading (or laboratory) axis (see Fig. 1) is accounted for through the rotation matrix Q in Eq. (12). For the 100 orientation Q is simply an identity. For the 110O1 orientation Q is   1 √ √1 0 2 2  − √1 √1 0  (36) 2 2 0 0 1 AN US For the 110O2 crystal orientation, the transformation matrix Q is  1  √ √1 0 2 2  0 0 1  1 1 √ − √2 0 2 For the 111O1 crystal orientation, the transformation matrixQ is  1  1 1 √ √ √ 3 √2 6 3 − √16   0 ED For the 111O2 crystal orientation, the transformation matrix Q is  1  1 1 CE PT √  − √31  2 √1 6 √ (38) √1 2 M  − √31  6 − √12 (37) 3 0 − √26 √ 3 √1 2 √1 6   (39) AC Table 2: Values of the material parameters for power law creep, Eq. (25), and the coefficient β characterizing the cubic plastic anisotropy, Eq. (23). η σ0 n f0 β2 1.53 × 10−9 sec−1 470.65 MPa 5 0.01 0.633 The values of material parameters, the initial void volume fraction and the anisotropy coefficients are tabulated in Table 2. The coefficient β can be expressed in terms of the relative creep resistance along two crystal orientations, for example any of the < 100 > and < 110 > directions, respectively. Alternatively, if experimental data are not available, β can 10 ACCEPTED MANUSCRIPT CR IP T be determined by adjusting its value such that the effective strain versus time response for two crystal orientations (under the same creep loading conditions) is in reasonable agreement with the corresponding response from finite element unit cell calculations. Here, we used the effective strain Ee versus time curves reported by Srivastava and Needleman (2015) for the 100, 110O1 and 111O1 crystal orientations for creep with χ = 3 and L = −1. Specifically, the values of η and n in Eq. (25) were taken the same as used in Srivastava and Needleman (2015) while σ0 and β were adjusted to obtain a best fit. Some calculations were repeated with other values of κ11 , κ12 , and κ44 and it was verified that the dependence on κij is only through the value of β 2 . Also, since the specific value of κ12 does not matter, identical values of β 2 , and therefore identical numerical results can be obtained with κ12 = 0. We note that specialization of the Hill (1948) criterion to cubic symmetry correspond to κ11 = 1 and κ12 = 0, so that the same results as presented here can be obtained using Benzerga and Besson (2001)2 potential specialized to cubic symmetry. PT ED M AN US 3.1. Constitutive Model Predictions (b) CE (a) AC Figure 2: Model predictions for the effect of crystal orientation on (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) for creep loading corresponding to three Lode parameter values, L = −1, 0 and 1 for a stress triaxiality value, χ = 3. In this subsection, we present the model predictions for the evolution of effective strain and porosity with time for various crystal orientations and creep loading conditions. The effective creep strain, Ee , is defined as r Z t 2 0 0 D D (40) Ee = Ėe (s)ds , Ėe = 3 ij ij 0 2 Note that Eqs. (19) and (54) of Benzerga and Besson (2001) coincide with Eqs. (24) and (22) if the values of their hi are taken as h1 = h2 = h3 = 1, h4 = h5 = h6 = β 2 . 11 ACCEPTED MANUSCRIPT 1 100, L=-1, 0, 1 110O1, L=1 0.8 3.5 110O1, L=-1 110O2, L=1 111O1, L=1 111O2, L=1 3 110O1, L=0 2.5 0.6 Ee 2 110O1, L=0 0.4 1.5 110O2, L= 0 111O1, L=-1 1 0 3 6 9 7 12 111O1, L=0 111O2, L=0 110O2, L= 0 111O1, L=-1 0.2 0 100, L=-1, 0, 1 110O1, L=1 110O1, L=-1 110O2, L=1 111O1, L=1 111O2, L=1 CR IP T Ee 111O1, L=0 111O2, L=0 15 0.5 18 Time (10 sec) 0 3 6 9 7 12 15 18 Time (10 sec) AN US (a) (b) Figure 3: Model predictions for the effect of crystal orientation on (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) for creep loading corresponding to three Lode parameter values, L = −1, 0 and 1 for a stress triaxiality value, χ = 2/3. AC CE PT ED M 0 is the deviatoric part of the rate of deformation tensor, Dij , given by Eq. (26). where Dij R Similarly, the void volume fraction, f , is f = f˙dt with f˙ given by Eq. (27). The predictions of the model for all five crystal orientations, Fig. 1, subjected to creep loading conditions corresponding to L = −1, 0 and 1 for χ = 3 are shown in Fig. 2. The material parameters used in these simulations are given in Table 2. First, the results for the 110O1 and 110O2 orientations with L = −1 must coincide, as is also the case for the 111O1 and 111O2 orientations, with L = −1 since the difference in these loadings corresponds to interchanging Σ02 and Σ03 ; and for loadings such L = −1, Σ02 and Σ03 are equal. Therefore, to limit the amount of text in the figures, the results for 110O2 and 111O2 with L = −1 are not marked. Also, as seen from Eqs. (A.3) and (A.17) the mechanical response for the 100 orientation (all values of L) and those for 110O1 with L = 1 is the same. Also, from Eqs. (A.15) and (A.19) follows that the mechanical responses for the 110O1 with L = −1 and 110O2 with L = 0 coincide. In addition, from Eqs. (A.22) and (A.25) the results for the 111O1 and 111O2 orientations coincide for both L = 0 and L = 1. Hence, the prediction is that for all three values of Lode parameter considered the response of the 111O1 and 111O2 orientations coincide, and this is seen in Figs. 2 and 3. In Fig. 2, with χ = 3, the 100 orientation for all values of Lode parameter L and the 110O1 orientation for L = 1 give the fastest growth rate for the creep strain Ee and the void volume fraction f . The 111O1 orientation for L=-1 and the 110O2 orientation for L = 0 give the slowest growth rates of Ee and f . The constitutive model also predicts that the evolution of Ee and f is the same for the 110O1 orientation with L = −1, the 110O2 orientation with L = 1, and both the 111O1 and 111O2 orientations with L = 1. 12 ACCEPTED MANUSCRIPT ED M AN US CR IP T Fig. 3 shows the predictions of the constitutive model for χ = 2/3. The ordering of the various orientations is the same but obviously the growth rates are very much reduced compared to the cases when the triaxiality is χ = 3. Calculations were also carried out with n = 7 but with all other material parameters as in Table 2 for all five crystal orientations considered, all three values of L and χ = 3. As expected, increasing n results in an increased creep rate and void growth rate as seen in Fig. 4. However, the dependence of creep deformation and porosity evolution on crystal orientation remains the same (i.e. the fastest creep rate and void growth rate corresponds to 100 and the slowest creep rate and void growth rate corresponds to 110O2 with L = 0). CE PT Figure 4: Model predictions for the effect of crystal orientation with a creep exponent n=7 on (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) for creep loading corresponding to three Lode parameter values, L = −1, 0 and 1 for a stress triaxiality value, χ = 3. AC The model predictions show the strong effect of crystal orientation on the creep response. Only for the 100 orientation, is the behavior of the single crystal is close to isotropic behavior in that the response for 100 is independent of the value of the Lode parameter L. The main difference between the response using isotropic Gurson potential and the single crystal potential for the 100 orientation is the value of the parameter h in the cosh term. For the Gurson potential h = 2, while for the single crystal characterized by the plastic coefficients given in Table 2, h = 2.3258. Thus, void growth is slower for the 100 single crystal as seen in Fig. 5b. Also, as seen in Fig. 5a, the evolution of creep strain is slower for the 100 single crystal. 13 CR IP T ACCEPTED MANUSCRIPT AN US (a) (b) M Figure 5: A comparison between the (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) for an isotropic material (β =1) and for a single crystal in the 100 orientation characterized by the same set of material parameters as in Table 2 when subject to creep loading corresponding to Lode parameter value, L = −1, for a stress triaxiality value, χ = 3. AC CE PT ED 3.2. Comparison with Cell Model Predictions Figs. 6 to 10 show the comparison of the constitutive model and cell model predictions. In viewing the comparison, it is important to note that the constitutive model only specifies the cubic symmetry. The constitutive relation used in the single crystal cell model calculations of Srivastava and Needleman (2015) specifies discrete slip systems. Furthermore, the cell model calculations were carried out within a finite deformation framework so that lattice rotations leading to slip system reorientation and void shape changes were accounted for. The cell model calculations also account for void-void interactions, which are not accounted for in the simple phenomenological constitutive model. Fig. 6 compares the constitutive model predictions with the unit cell model results of Srivastava and Needleman (2015) for the 100 orientation with L = −1, 0, 1 and the 110O1 orientation with L = 1 for a triaxiality value χ = 3. Consistent with the constitutive model predictions the cell model response is identical for these orientations. Also, up to an overall strain of ≈ 0.05, there is very good quantitative agreement between the constitutive model and cell model predictions for the evolution of overall creep strain Ee in Fig. 6a. The constitutive model prediction for the porosity evolution in Fig. 6b underestimates the cell model predictions except for very small amounts of void growth. Creep curves and porosity evolution for the 110O1 orientation with L = −1, the 110O2 orientation with L = 1, the 111O1 orientation with L = 1 and the 111O2 orientation with L = 1 with χ = 3 are shown in Fig. 7. As predicted by the constitutive model, the cell 14 CR IP T ACCEPTED MANUSCRIPT AN US (a) (b) M Figure 6: Comparison between the (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) as predicted by the model (line) and the unit cell calculations of Srivastava and Needleman (2015) (symbols) for 100 crystal orientation for creep loading corresponding to L = −1, 0 and 1 for χ = 3, and for 110O1 crystal orientation for creep loading corresponding to L = 1 for χ = 3. AC CE PT ED model responses for these orientations and Lode parameter values are essentially identical at small strains. Furthermore, for these orientations there is good quantitative agreement between the constitutive model and cell model predictions for both the evolution of creep strain Ee and porosity f in the early stages of deformation. At larger strains, the cell model predictions for the responses of these orientations separates mainly due to an evolving slip mode due to the reorientation of slip systems with increasing deformation. The predictions of the constitutive model are compared with those of the cell model calculations for the 111O1 and 111O2 orientations with L = 0 in Fig. 8. As in Fig. 7, the quantitative as well as qualitative agreement is very good at small strains and the qualitative trends are well represented by the constitutive model. Fig. 9 compares the constitutive model predictions with the cell model predictions for the 110O2 orientation with L = 0 and the 111O1 orientation with L = −1. The constitutive model predicts that the responses in these two cases are identical. Fig. 9a shows that the cell model predictions for the evolution of the creep strain Ee for these two cases deviates at small strains. However, the porosity evolution for these two cases is essentially identical for relatively large amounts of void growth, up to f /f0 ≈ 4, and there is very good quantitative agreement between the constitutive model and cell model predictions over this range. A comparison of constitutive model and cell model predictions with χ = 2/3 is shown in Fig. 10. The cases considered are the 100, 110O1 and 111O1 orientations with L = −1 for all three orientations. The quantitative agreement is not as good as for χ = 3 but the 15 CR IP T ACCEPTED MANUSCRIPT AN US (a) (b) AC CE PT ED M Figure 7: Comparison between the (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) as predicted by the model (line) and the unit cell calculations of Srivastava and Needleman (2015) (symbols) for 110O1, 110O2, 111O1 and 111O2 crystal orientations for creep loading corresponding to L = −1, 1, 1 and 1, respectively, for χ = 3. (a) (b) Figure 8: Comparison between the (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) as predicted by the model (line) and the unit cell calculations of Srivastava and Needleman (2015) (symbols) for 111O1 and 111O2 crystal orientations for creep loading corresponding to L = 0 for χ = 3. 16 CR IP T ACCEPTED MANUSCRIPT AN US (a) (b) AC CE PT ED M Figure 9: Comparison between the (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) as predicted by the model (line) and the unit cell calculations of Srivastava and Needleman (2015) (symbols) for 110O2 and 111O1 crystal orientations for creep loading corresponding to L = 0 and 1, respectively, for χ = 3. (a) (b) Figure 10: Comparison between the (a) creep curves (effective creep strain, Ee , versus time) and (b) porosity evolution (relative void volume fraction, f /f0 , versus time) as predicted by the model (line) and the unit cell calculations of Srivastava and Needleman (2015) (symbols) for 100, 110O1 and 111O1 crystal orientations for creep loading corresponding to L = −1 for χ = 2/3. 17 ACCEPTED MANUSCRIPT qualitative trends seen in the numerical results are mostly consistent with the constitutive model predictions except for the relative growth of the void volume fraction for the 110O1 and 111O1 orientations where the cell model and constitutive model results differ regarding the orientation with the faster growth rate of porosity. CR IP T 4. Discussion Table 3: Values of the Sachs factor, Sxtal , reproduced from Srivastava and Needleman (2015) and the values of the Sachs factor, Stheory estimated using the phenomenological theory for all the orientations and stress states analyzed. Stheory 10.28 10.28 10.28 3.91 7.70 10.28 2.60 3.91 2.60 3.00 3.91 3.00 3.91 ED M AN US Orientation L Sxtal 100 −1 9.98 100 0 12.2 100 +1 9.98 110O1 −1 4.99 110O1 0 12.0 110O1 +1 9.98 110O2 0 0.37 110O2 +1 4.99 111O1 −1 0.66 111O1 0 2.09 111O1 +1 5.02 111O2 0 1.13 111O2 +1 2.98 AC CE PT In Srivastava and Needleman (2015) the qualitative effect of crystal orientation on creep curve for various values of the Lode parameter was shown to be revealed by the Sachs factor. The Sachs factor is essentially a normalized value of the plastic dissipation accompanying creep. The Sachs factor S is given by Σij Dij S= (41) Ẇref where Ẇref is some convenient normalizing quantity. In Srivastava and Needleman (2015) Ẇref was taken to be γ̇0 Σe , with Σe the imposed Mises equivalent stress and γ̇0 a reference strain rate. Similarly for the phenomenological theory Ẇref is taken as ηΣe , with Σe the imposed Mises equivalent stress and η the reference creep rate in Eq. (25). For the fully dense rigid plastic matrix material the Sachs parameter is independent of the hydrostatic stress. The Sachs factor, Stheory , estimated using the phenomenological theory for the fully dense single-crystal (f0 = 0) for the orientations and stress states analyzed together with the Sachs factor, Sxtal from Srivastava and Needleman (2015) are tabulated in Table 3. The larger 18 ACCEPTED MANUSCRIPT CE PT ED M AN US CR IP T the value of S the softer the creep response, i.e. the greater the creep rate. Comparing the constitutive model predictions in Figs. 2a and 3a with the values of S in Table 3 shows that the constitutive model predictions correlate well with the qualitative difference in the values of the Sachs factor. For example, in Fig. 2a the creep rate for the 100 orientation for all three values of L are roughly the same, consistent with the values of Stheory equal to 10.28 or Sxtal varying from 9.98 to 12.2. For 110O1 orientation the creep rates for L = 0 and L = 1 are greater than for L = −1 as predicted by the model and the Sachs factor, while with the 110O2 orientation the creep rate is largest for L = ±1 and lowest for L = 0. Also for the 111O1 and 111O2 orientations, the ordering of the creep responses predicted by the model is, to an extent, consistent with the values of the Sachs factors in Table 3. While the flow potential for the fully-dense single-crystal (see Eqs. 13 and 24 with f ≡ 0) does not account for the geometry of the slip systems, it satisfies the requirements of form-invariance with respect to the transformations belonging to the group of symmetries of the fcc single-crystal i.e. it accounts for the fact that the single crystal exhibits three orthogonal mirror planes [100], a three-fold rotation axis in the < 111 > direction, and a mirror plane at 45o with respect to [100]. Thus, the flow potential, and hence the evolution of porosity and of plastic deformation, account for the influence of the basic single-crystal orientation with respect to the loading. For a crystal there are two sources of anisotropy, one is the basic crystal structure, cubic for the crystal structure considered, and the other is the slip system anisotropy. The flow potential for the fully dense solid, Eq. (13) with f = 0, accounts for the former but not for the latter. Our results in Figs. 6 to 9 for a stress triaxiality value χ = 3 show that, in the early stages of deformation, there is very good agreement between the phenomenological theory and crystal plasticity finite element calculations with unit-cells containing a single void (with the exception of the creep curves for 111 and L = 0). Also, the plastic potential, Eq. (13), reduces to the original Gurson (1975) potential for an isotropic matrix where it is known that the q parameters introduced by Tvergaard (1981, 1982) are needed to obtain good agreement with cell model calculations. The divergence of the phenomenological predictions from the crystal plasticity cell model predictions at later times is most likely due to the evolution of slip system activity and void shape changes. The much less good agreement for the lower values of stress triaxiality in Fig. 10 is expected for any Gurson type relation that neglects void shape changes. AC 5. Conclusions The predictions of the proposed phenomenological model for the effect of crystal orientation and applied stress state (characterized by the value of stress triaxiality and Lode parameter) on the evolution of creep strain strain and porosity are in good qualitative agreement with the finite strain unit cell model calculations for fcc single crystals. Consistent with the cell model calculations, for the 100 crystal orientation that gives rise to nearly isotropic response, the phenomenological model predicts no effect of the Lode parameter on the dilatational response while for anisotropic crystal orientations a significant influence of the Lode parameter is predicted. 19 ACCEPTED MANUSCRIPT CR IP T While the very strong influence of crystal anisotropy and loading path on the evolution of the porosity and plastic deformation has been seen in previous numerical finite element cell model studies, only partial qualitative explanations of the numerical results have been presented (Yerra et al., 2010; Srivastava and Needleman, 2015). Our model predictions have revealed previously unrecognized features of the creep response of porous single crystals. For example, results obtained from the model reveal that the creep response is the same for certain crystal orientations and loadings. Our results show that: • The creep response for the 111 crystal orientation under axisymmetric loading corresponding to a Lode parameter value of L = −1 should be the same as that for 110O2 crystal orientation under loading corresponding to L = 0. AN US • The creep response for the 111 crystal orientation under axisymmetric loading corresponding to L = 1, should be the same as that of the 110 crystal orientation under axisymmetric loading corresponding to L = −1, and that of the 110O2 crystal orientation under axisymmetric loading corresponding to L=1; • The creep response for the 110O1 crystal orientation under axisymmetric loading corresponding to L = 1 should be the same as that of the [100] crystal orientation for all loadings. AC CE PT ED M Our model calculations also predict that the slowest void growth rate should correspond to loading with a Lode parameter value L = 0 for the 110O2 orientation (and to L = −1 for the 111 orientation), while the fastest rate of void growth should correspond to loading with L = 1 for the 110O1 orientation (and for the 100 orientation for all values of L). On the other hand, the rate of void growth for a loading with L = 0 in the 111 orientation is predicted to be less than the rate of void growth for axisymmetric loadings in the 110 orientation and in the 100 orientation. We have applied the analytical model in Eq. (1) to the description of creep of porous fcc crystals that deform by crystallographic slip described by Schmid law, so that the parameter k in Eq. (1) was set to zero. However, it is known, see for example Vitek et al. (2004), that for certain bcc single crystals (e.g. molybdenum), there is a tension-compression asymmetry due to non-planar spreading of individual dislocations. For such a bcc porous single-crystal a value of k different from zero should be used. Appendix In this Appendix we derive simple explicit expressions for the effective stress σ̃e2 for all orientations and values of Lode parameter presented in the numerical results. These expressions follow from the fact that for cubic symmetry there is a simple expression for the effective stress, Eq. (24), valid for all orientations and Lode parameter values that only depends on the parameter β 2 , defined in Eq. (23), and not on the individual values of κij . 20 ACCEPTED MANUSCRIPT 0 0 0 0 0 0 From Eqs. (16) to (19), σ̃12 = σ̃21 , σ̃13 = σ̃31 , σ̃23 = σ̃32 , using the definition of σ̃e2 0 0 0 Eq. (14), and σ11 + σ22 + σ33 = 0 gives 

 02 02 02 02 02 02 + σ13 + σ23 + σ11 + σ22 + σ33 (κ11 − κ12 )2 (A.1) σ̃e2 = m̂2 2κ244 σ12 Substituting the definition of m̂2 , Eq. (21), into Eq. (A.1) gives Eq. (24). AN US CR IP T Effect of Orientation The deviatoric components of the Cauchy stress tensor on the crystal axes, σij0 , are related to the applied (along laboratory axes) stress deviator, Σ0ij through Eq. (12). Hence for the orientations considered here, see Fig. 1, explicit expressions can be given for σ̃e2 in terms of the applied stress deviator components Σ01 , Σ02 and Σ03 . For the 100 orientation  0  Σ1 0 0 σ 0 =  0 Σ02 0  (A.2) 0 0 0 Σ3 and  3  02 02 2 Σ1 + Σ02 (A.3) 2 + Σ3 = Σe 2 Hence, for the 001 orientation the plastic potential φ and its derivatives are independent of the value of the Lode parameter. For the 110O1 orientation   1 0 (Σ1 + Σ02 ) 12 (Σ01 − Σ02 ) 0 2    1 0  0 1 0 0 0  (Σ − Σ ) (Σ + Σ ) 0 (A.4) σ = 1 2 1 2 2 2     0 0 Σ03 PT ED M σ̃e2 = and CE σ̃e2 AC For the 110O2 orientation and    σ =   0 σ̃e2   3 3 02 β 2 0 0 2 = Σ + (Σ1 − Σ2 ) 2 2 3 2 1 2 (Σ01 + Σ03 ) 1 2 (Σ01 − Σ03 ) 1 2 (Σ01 1 2 (Σ01 − Σ03 ) 0 + 0 Σ03 ) 0    0    Σ02   3 3 02 β 2 0 0 2 = Σ + (Σ1 − Σ3 ) 2 2 2 2 21 (A.5) (A.6) (A.7) ACCEPTED MANUSCRIPT 1 " Σ01 2Σ02 + 3 3 2 +2 Σ01 3  1 3 σ̃e2 = 2 " Σ01 2Σ03 + 3 3 2 +2 3 Σ01 3  Σ02 6 Σ01 3 Σ02 3 − Σ01 3 Σ01 3 Σ03 2 + Σ01 3 − Σ02 2 − Σ02 2 + Σ03 6 Σ03 3 + 2 Σ01 3 Σ01 3 Σ03 6 Σ01 3 Σ01 Σ02 Σ03 + + 3 2 6 ED and + − Σ03 2 Σ01 Σ02 Σ03 + + 3 6 2 For the 111O2 orientation  Σ0     σ0 =     Σ01 3 + − + − Σ02 3 Σ01 3 2 Σ02 3 Σ02 3 Σ01 3 Σ01 3 + 2β 2 2 − Σ03 3 + Σ01 3 +  Σ02 6 − Σ02 6 − Σ03 2 Σ02 3 + Σ03 2 Σ01 Σ02 − 3 3 − Σ02 2 + Σ03 6 AN US 3 σ̃e2 = 2 3 Σ02 6 M and     σ0 =     + 2 + − 2Σ03 3 Σ03 3 Σ01 3 Σ01 3 + 2β 2 2 +           (A.8) CR IP T For the 111O1 orientation  Σ0 − Σ02 2 Σ03 3 + Σ03 6 Σ01 Σ03 − 3 3 2 +  2 !# Σ01 Σ02 Σ03 + − 3 6 2 (A.9)          2 (A.10) +  2 !# Σ01 Σ02 Σ03 − + 3 2 6 (A.11) AC CE PT Effect of the Lode Parameter For L = −1 the applied stress deviator the principle values of the stress deviator Σ0 are Σ01 , −Σ01 /2, −Σ01 /2; for L=0 these values are Σ01 , 0, −Σ01 ; and for L=1 the principle values are Σ01 , Σ01 , −2Σ01 . It is convenient to write the components of Σ0 in terms of the applied effective stress Σe especially in the present context since the calculations are carried out for cases for which the value of Σe , is the same for each value of the Lode parameter L. The matrix of components of Σ0 for each value of the Lode parameter L are then given by L = −1,  2  Σ 0 0 3 e 0  Σ0 =  0 − 13 Σe (A.12) 1 0 0 − 3 Σe L = 0,  Σ0 =  √1 Σe 3 0 0 22  0 0  0 0 0 − √13 Σe (A.13) ACCEPTED MANUSCRIPT L = 1,   0 0 1 Σ 0  3 e 2 0 − 3 Σe 1 Σ 3 e Σ0 =  0 0 (A.14) For the 110O2 orientation: σ̃e2 = Σ2e  for L = 1 AN US σ̃e2 = Σ2e CR IP T Variation of σ̃e2 with orientation and Lode parameter value For the 100 orientation, σ̃e2 = Σ2e for all values of L as seen in Eq. (A.3). For the 110O1 orientation:   3 2 2 2 1 σ̃e = Σe + β for L = −1 4 4   1 2 2 3 2 for L = 0 σ̃e = Σe + β 4 4 1 3 2 + β 4 4  for L = −1 ED For the 111O1 orientation: M σ̃e2 = β 2 Σ2e for L = 0   3 2 2 2 1 + β σ̃e = Σe for L = 1 4 4 CE PT σ̃e2 = β 2 Σ2e for L = −1   11 2 1 2 2 + β σ̃e = Σe for L = 0 12 12   3 2 2 2 1 σ̃e = Σe + β for L = 1 4 4 (A.15) (A.16) (A.17) (A.18) (A.19) (A.20) (A.21) (A.22) (A.23) AC For the 111O2 orientation: σ̃e2 = β 2 Σ2e for L = −1   11 2 1 2 2 σ̃e = Σe + β for L = 0 12 12   3 2 2 2 1 σ̃e = Σe + β for L = 1 4 4 (A.24) (A.25) (A.26) Note that in all cases σ̃e2 = Σ2e for β 2 = 1 (which corresponds to an isotropic solid). 23 ACCEPTED MANUSCRIPT References AC CE PT ED M AN US CR IP T Ahzi, S., Schoenfeld, S.E., 1998. Mechanics of porous polycrystals: a fully anisotropic flow potential. International Journal of Plasticity, 14, 829-839. Arminjon, M., 1991. A regular form of the Schmid law. Application to the ambiguity problem. Textures and Microstructures, 14, 1121-1128. Asaro, R.J., Needleman, A., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metallurgica 33, 923-953. Benzerga, A.A., Besson, J., 2001. Plastic potentials for anisotropic porous solids. European Journal of Mechanics-A/Solids, 20, 397-434. Benzerga, A., Leblond, J-B., 2010. Ductile fracture by void growth to coalescence. Advances in Applied Mechanics, 44, 169-305. Bishop, J. F. W., Hill, R., 1951. A theory of the plastic distortion of a polycrystalline aggregate under combined stresses. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 42, 414-427. Cazacu, O., Plunkett, B., Barlat, F., 2006. Orthotropic yield criterion for hexagonal closed packed metals. International Journal of Plasticity, 22, 1171-1194. Cazacu, O., Ionescu, I. R., Yoon, J. W., 2010. Orthotropic strain rate potential for the description of anisotropy in tension and compression of metals. International Journal of Plasticity, 26, 887-904. Cazacu, O., Stewart, J. B., 2013. Analytical criterion for porous solids containing cylindrical voids in an incompressible matrix exhibiting tensioncompression asymmetry. Philosophical Magazine, 93, 1520-1548. Gurson, A.L., 1975. Plastic flow and fracture behavior of ductile materials incorporating void nucleation, growth and interaction. Ph.D. thesis, Brown University, Providence, RI. Ha, S., Kim, K., 2010. Void growth and coalescence in fcc single crystals. International Journal of Mechanical Sciences, 52, 863-873. Han, X., Besson, J., Forest, S., Tanguy, B., Bugat, S., 2013. A yield function for single crystal containing voids. International Journal of Solids and Structures, 50, 2115-2131. Hill, R., 1948. A theory of yielding and plastic flow of anisotropic solids. Proceedings of the Royal Society of London A 193, 281-297. Hill, R., 1967. The essential structure of constitutive laws for metal composites and polycrystals. Journal of the Mechanics and Physics of Solids, 15, 79-95. Idiart, M. I., Castaeda, P. P., 2007. Variational linear comparison bounds for nonlinear composites with anisotropic phases. II. Crystalline materials. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 463, 925-943. Lebensohn, R.A., Cazacu, O., 2012. Effect of single-crystal plastic deformation mechanisms on the dilatational plastic response of porous polycrystals. International Journal of Solids and Structures, 49, 3838-3852. Mandel, J., 1972. Plasticité classique et viscoplasticité. CISM Courses and Lectures, Vol.97, International Center for Mechanical Sciences, Springer-Verlag, Wien-New York. Mbiakop, A., Constantinescu, A., Danas, K., 2015. A model for porous single crystals with cylindrical voids of elliptical cross-section. International Journal of Solids and Structures, 64, 100-119. Monchiet, V., Cazacu, O., Charkaluk, E., Kondo, D., 2008. Macroscopic yield criteria for plastic anisotropic materials containing spheroidal voids. International Journal of Plasticity, 24, 1158-1189. Morin, L., Madou, K., Leblond, J.B., Kondo, D., 2014. A new technique for finite element limit-analysis of Hill materials, with an application to the assessment of criteria for anisotropic plastic porous solids. International Journal of Engineering Science, 74, 65-79. Needleman, A., Tvergaard, V., Hutchinson, J.W., 1992. Void growth in plastic solids. Topics in Fracture and Fatigue, (ed. by A.S. Argon), Springer-Verlag, New York, 145-178. Pan, J., Saje, M., Needleman, A., 1983. Localization of deformation in rate sensitive porous plastic solids. International Journal of Fracture, 21, 261-278. Paux, J., Morin, L., Brenner, R., Kondo, D., 2014. An approximate yield criterion for porous single crystals. European Journal of Mechanics-A/Solids, 51, 1-10. 24 ACCEPTED MANUSCRIPT AC CE PT ED M AN US CR IP T Ponte Castaneda, P., 2002. Second-order homogenization estimates for nonlinear composites incorporating field fluctuations: Itheory. Journal of the Mechanics and Physics of Solids, 50, 737-757. Rice, J.R., 1970, On the structure of stress-strain relations for time-dependent plastic deformation in metals. Journal of Applied Mechanics, 37, 728-737. Srivastava, A., Needleman, A., 2012. Porosity evolution in a creeping single crystal. Modelling and Simulation in Materials Science and Engineering, 20, 035010. Srivastava, A., Needleman, A., 2013. Void growth versus void collapse in a creeping single crystal. Journal of the Mechanics and Physics of Solids, 61, 1169-1184. Srivastava, A. and Needleman, A., 2015. Effect of crystal orientation on porosity evolution in a creeping single crystal. Mechanics of Materials, 90, 10-29. Stewart, J. B., Cazacu, O., 2011. Analytical yield criterion for an anisotropic material containing spherical voids and exhibiting tension-compression asymmetry. International Journal of Solids and Structures, 48, 357-373. Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. International Journal of Fracture, 17, 389-407. Tvergaard, V., 1982. On localization in ductile materials containing spherical voids. International Journal of Fracture, 18, 237-252. Tvergaard, V., 1990. Material failure by void growth to coalescence. Advances in Applied Mechanics, 27, 83-151. Wan, J.S., Yue, Z.F., Lu, Z.Z., 2005. Casting microporosity growth in single-crystal superalloys by a three dimensional unit cell analysis. Modelling and Simulation in Material Science and Engineering, 13, 875-892. Yerra, S.K., Tekoglu, C., Scheyvaerts, F., Delannay, L., Van Houtte, P., Pardoen, T., 2010. Void growth and coalescence in single crystals, International Journal of Solids and Structures. 47, 1016-1029. Yu, Q.M., Hou, N.X., Yue, Z.F., 2010. Finite element analysis of void growth behavior in nickle-based single crystal superalloys, Comutational Materials Science. 48, 597-608. Vitek, V., Mrovec, M., Bassani, J. L., 2004. Influence of non-glide stresses on plastic flow: from atomistic to continuum modeling. Materials Science and Engineering: A, 365, 31-37. Worswick, M.J., Wong, B., Pick, R.J., 1991. Void growth during high velocity impact: experiment and model. Le Journal de Physique IV, 1(C3), C3-605. 25