Finite element simulation of ductile fracture in polycrystalline materials using a regularized porous crystal plasticity model

In the present study, a hypoelastic–plastic formulation of porous crystal plasticity with a regularized version of Schmid’s law is proposed. The equation describing the effect of the voids on plasticity is modified to allow for an explicit analytical solution for the effective resolved shear stress. The regularized porous crystal plasticity model is implemented as a material model in a finite element code using the cutting plane algorithm. Fracture is described by element erosion at a critical porosity. The proposed model is used for two test cases of two- and three-dimensional polycrystals deformed in tension until full fracture is achieved. The simulations demonstrate the capability of the proposed model to account for the interaction between different modes of strain localization, such as shear bands and necking, and the initiation and propagation of ductile fracture in large scale polycrystal models with detailed grain description and realistic boundary conditions.


Introduction
For many types of Al alloys and steels fracture occurs by ductile fracture mechanisms, namely the nucleation, growth and coalescence of microscopic voids. Ductile fracture has been studied extensively due to its practical importance, but also because it is a fundamentally interesting complex multiscale process. The ductile fracture process involves material deforming plastically around micron-sized particles and voids at complex stress states.
The problem of plastic deformation in a material containing a void has been approached analytically since McClintock's analysis of a cylindrical void in an infinite plastic medium (McClintock 1968) and Rice and Tracey's expression for exponential void growth due to the hydrostatic stress (Rice and Tracey 1969). In the seminal work of Gurson (1977), limit analysis was applied to a spherical void in a von Mises plastic medium to derive the yield function for the porous material. A phenomenological approach based on continuum thermodynamics was proposed in Rousselier (1981), where ductile damage was introduced as a pressure-dependent term in the macroscopic plastic potential. At the same time, the development of the finite element method allowed for modelling the behaviour of a void in a medium numerically. One of the first finite element models of a porous material unit cell was presented by Needleman (1972). The unit cell models may predict the behaviour of any complex non-linear material with voids of arbitrary shape, size and position relative to each other. It can also provide both global averaged stresses and strains, and the local microscopic stress and strain fields around the voids, including the effects of void interaction and strain localization (Daehli et al. 2017a, b). Unit cell simulations provide a way to validate the analytical porous plasticity models (Guo and Li 2019). The Gurson model showed some deviations from the unit cell results and accordingly a phenomenological correction to the model was proposed in Tvergaard (1981). Coefficients were introduced to account for the void shape change from spherical during deformation, void interactions and work hardening of the material. Other modifications to the Gurson model were proposed over the years. Several populations of voids with different sizes were analysed in Perrin and Leblond (1990) and void coalescence in Koplik and Needleman (1988). Other modifications include improved formulation for matrix material with high work hardening rate (Leblond et al. 1995), kinematic hardening (Besson and Guillemer-Neel 2003;Mear and Hutchinson 1985;Mühlich and Brocks 2003), ellipsoidal void shape (Benzerga et al. 2004;Gologanu et al. 1993;Pardoen and Hutchinson 2000), and viscoplastic matrix (Flandi and Leblond 2005;Moran et al. 1991).
Thermomechanical processing of metals may reorient the grains and produce certain crystallographic textures. Texture is the primary source of plastic anisotropy in metals (Barlat 1987). Several works combined the porous plasticity model of Gurson with plastic anisotropy of the matrix material. One simple way to obtain such a modification is to replace the von Mises norm of the stress tensor in the Gurson yield function with the norm of the corresponding anisotropic yield function. Most often the anisotropic yield function by Hill (1948) has been used (Doege et al. 1995;Grange et al. 2000;Rivalin et al. 2001;Wang et al. 2004). The anisotropic yield function by Barlat et al. (2005) was used in Daehli et al. (2017a) instead. On the other hand in Benzerga and Besson (2001), Benzerga et al. (2004) and Monchiet et al. (2008) the porous plasticity model is derived from the initial assumption of an anisotropic matrix.
A more direct way to include plastic anisotropy into the modelling of ductile materials is using crystal plasticity instead of phenomenological anisotropic plasticity. In the case of ductile fracture in metals, the voids nucleate and grow inside grains or at the boundaries between two grains at the micron and submicron scale. Therefore, it is natural to try using crystal plasticity to model the ductile crystal grain behaviour instead of phenomenological plasticity, which describes the averaged plastic behaviour of polycrystals. An early attempt to analyse the behaviour of voids inside a crystal grain can be found in Nemat-Nasser and Hori (1987) for a 2D case. Analytical derivation of a porous plasticity model is not trivial even for a von Mises material. It is even more challenging for the complex plasticity description of crystals. Several approaches exist that allowed such derivations. One approach uses limit analysis, in the same vein as Gurson (1977). An example of this approach is presented in Han et al. (2013), which also adds the additional free parameters to fit the behaviour of the model to unit cell simulations, analogous to Tvergaard's extension of the Gurson model. The model was reformulated for finite deformations by Ling et al. (2016). Other examples of this approach include Paux et al. (2015) and Paux et al. (2018). Another approach, called the non-linear variational homogenization method, was developed by Ponte Castañeda (1991). It is a general approach of modelling heterogeneous plastic materials, which was subsequently applied to voided crystals with ellipsoidal voids (Mbiakop et al. 2015a, b). Interestingly, when this approach is combined with unit cell simulation results, it produces yield functions that are quite similar to the one derived in Han et al. (2013). In recent works (Song and Ponte Castañeda 2017a, b, c), the evolution of the void shape and the subsequent morphological anisotropy of the porous material were included in the model. A third approach, based on sequential laminates of infinite rank, was used to derive the yield function of a porous crystal in Joëssel et al. (2018). Another approach to studying fracture in polycrystals, which is currently developing, is the combination of phase field modelling and crystal plasticity (De Lorenzis et al. 2016;Diehl et al. 2017;Padilla and Markert 2017). In addition, some works try to formulate and implement micromorphic size-dependent porous crystal plasticity (Ling et al. 2018) or add a void coalescence criterion (Hure 2019;Siddiq 2019).
In these studies, the porous crystal plasticity models were developed and, in some cases, validated using unit cell models. Very few works have attempted to further use the models to investigate fracture phenomena in crystals, although in Ling et al. (2016) and Ling et al. (2018) 2D single crystal tensile tests were simulated. The crystal plasticity models are computationally heavy by themselves, and introducing porous plasticity introduces another layer of complexity. Porous crystal plasticity models derived based on limit analysis (Han et al. 2013) or variational homogenization (Song and Ponte Castañeda 2017a) require iterative calculations in each timestep of the temporal integration of the rate constitutive equations. On the other hand fracture in a polycrystal is preceded by a complex stress and strain history with a combination of different types of strain localization (necking and shear banding), crystal orientation evolution and heterogeneous stress fields (Di Gioacchino and da Fonseca 2015;Guery et al. 2016;Lim et al. 2014). Modelling these processes requires finite element models of the grain structure with a high resolution that are stable for large deformation processes. Therefore, an accurate, robust and computationally effective numerical implementation of porous crystal plasticity models for finite element analysis is an important task.
In the present study, an implementation of the porous crystal plasticity model in the finite element solver LS-DYNA is described, where the regularization of Schmid's law proposed by Zamiri and Pourboghrat (2010) is adopted for the porous crystal plasticity model proposed by Han et al. (2013). The main equation of the porous crystal plasticity model is modified to allow explicit solution for the effective resolved shear stress. The proposed constitutive model is implemented in the commercial finite element code with an element erosion criterion based on critical porosity. It is then used to simulate the deformation until fracture of polycrystals in plane-strain tension (2D) and uniaxial tension (3D), demonstrating the complex interaction between the heterogeneous grain structure, different modes of strain localization and fracture initiation and propagation in polycrystals with realistic grain structures and boundary conditions.

Regularized porous crystal plasticity model
The plastic deformation of crystals is assumed to be due to plastic slip on a set of crystallographic slip systems, defined by the slip plane normal n a ð Þ and slip direction m a ð Þ , where n a ð Þ and m a ð Þ are unit vectors, a 2 1; . . .; N ½ signify the slip system, and N is the number of slip systems. Twinning, grain boundary sliding, and other deformation types are not considered. The FCC lattice is assumed, with N ¼ 12 independent slip systems of the 111 f gh110i family. The elastic deformations are assumed to remain small, while the plastic deformations and rotations can be finite. For a detailed description of crystal plasticity, the reader is referred to Roters et al. (2010).
The velocity gradient tensor L of the crystalline material can be additively decomposed into elastic L e and plastic L p parts: where F is the deformation gradient tensor. The velocity gradient tensor may be decomposed into the symmetric rate-of-deformation tensor D and the skewsymmetric spin tensor W: The and operations produce correspondingly the symmetric and skew-symmetric parts of the tensor. The total rate-of-deformation and spin tensors can also be decomposed into elastic and plastic parts: It will be assumed here that the elastic deformations are infinitesimal, which is a reasonable assumption for metals. Accordingly, the elastic spin tensor W e consists of an infinitesimal elastic contribution and rigid spin of the crystal lattice, whereas the plastic spin tensor W p is caused by plastic slip.
The crystal lattice undergoes finite rotations during the deformation. Therefore, it is convenient to define a co-rotational coordinate system that rotates with the crystal lattice. The orthogonal rotation tensor R defines the rotation from the global (fixed) basis e i to the co-rotational basisê i . In the following, denotes vectors and tensors defined in the co-rotational system: where R T Á R ¼ R Á R T ¼ I and I is the second-order identity tensor. By definition, the plastic spin W p does not affect the slip system vectors n a ð Þ and m a ð Þ , which rotate with the elastic spin W e : Accordingly, the rotation tensor evolves with the elastic spin: The slip system vectors in the global system are then related to the same vectors in the co-rotational coordinate system as: Similarly, the rate-of-deformation and spin tensors in the two basis systems are related as: The inverse transformations from the global to the co-rotational system are readily obtained by the orthogonality of the rotation tensor. A description of the co-rotational formulation of hypoelastic crystal plasticity may be found in Zhang et al. (2014).
The co-rotational stress tensorr is obtained from the rate form of the generalized Hooke's law in the corotational coordinate system: whereĈ is the fourth-order elasticity tensor in the corotational system, which may be assumed constant, andD e is the elastic part ofD. For the FCC lattice, the elasticity tensor has three independent constants,ĉ 11 , c 12 andĉ 44 , which describe the elastic anisotropy of the crystal. In Voigt notationĈ may be written asĈ The resolved shear stress s a is defined as the projection of the co-rotational stress tensorr onto slip system a: The porous crystal plasticity model proposed by Han et al. (2013) and reformulated for finite deformations in Ling et al. (2016) is applied in this study, using the co-rotational stress formulation. The void volume fraction denoted f is the evolving material parameter. The effective resolved shear stress s a eff on slip system a, which accounts for the effects of the voids, is defined by: where a, q 1 and q 2 are the parameters, like those introduced into the Gurson model by Tvergaard (1981), that improve the global accuracy and bring the predictions closer to unit cell simulation results. The von Mises norm of the Cauchy stress is denoted r vM and the hydrostatic stress is denoted r H . The sign of the effective resolved shear stress s a eff is the same as the sign of the corresponding resolved shear stress s a on slip system a.
An iterative process, e.g. a Newton-Raphson scheme, is required to find s a eff for each time step because Eq. (12) cannot be solved analytically, and this is considered a disadvantage in the numerical implementation of the porous crystal plasticity model. To circumvent this problem, the hyperbolic cosine in Eq. (12) is here approximated by the first four members of the Taylor polynomial: The polynomial approximation allows solving for s a eff explicitly. The resulting expression is a quartic equation with respect to 1=s a eff À Á 2 . Using more terms in the series expansion leads to higher-order polynomial equations, which have no known analytical roots.
Using less terms would result in less computations but becomes a compromise with accuracy. In Zamiri and Pourboghrat (2010), a regularized yield function for crystal plasticity is proposed, which is used in this study. The problem of elastoplastic deformation is interpreted as a constrained optimization problem, where the slip increments (or slip rates) are the Lagrange multipliers and the constraints are defined by: According to Kreisselmeier and Steinhauser (1980), the set of the 12 exact constraints can be replaced by the following form also known as the KSfunction: which is a domain enclosed by a smooth envelope of the convex polytope. The parameter q defines how close the envelope approaches the polytope. The plastic velocity gradient tensorL p is obtained from the yield function and the normality rule: where _ k is the plastic parameter and the factor 1 À f reflects the volume fraction of the single crystal in which plastic dissipation occurs (Besson 2009). Note that the symmetry ofr was not enforced in the differentiation to obtain the generally non-symmetric tensorL p by means of the normality rule, see Ling et al. (2016). The plastic rate-of-deformation tensorD p and the spin tensorŴ p are then defined as the symmetric and skew-symmetric parts ofL p , respectively. Using the chain rule, we get:L where the slip rates _ c a are expressed as: The partial derivative os a eff =or for the adopted polynomial approximation is given in Appendix 1. The plastic dissipation of the porous single crystal takes the form: where it was used that s a eff is a homogeneous function of degree one inr (Han et al. 2013). We note that the plastic dissipation is non-negative for f 1 and s a cr [ 0. By assuming a plastically incompressible matrix material, the evolution of the void volume fraction due to growth of existing voids is given by: In the current version of the regularized porous crystal plasticity model, we have neglected void nucleation and void coalescence, but these phenomena can be readily included at the cost of adding some extra material parameters.
The crystals work-harden during the plastic deformation due to dislocation accumulation. This is reflected in the model by the slip resistance increasing with the accumulated slip. The evolution of the slip resistance is described by the equation: where h ab is the instantaneous hardening matrix.
Various work hardening rules for single crystals are described in the literature. In the present study, the exact form of the work hardening rule is not important, thus for simplicity the Voce hardening rule is used.
The instantaneous hardening matrix is decomposed into the latent hardening matrix q ab and the selfhardening rate of slip systems H C ð Þ: where C is the total accumulated slip defined as: The latent hardening matrix is defined by: where a commonly used value of q ¼ 1:4 is assumed.
The self-hardening rate is equal to: where N V is the number of terms describing the work hardening, h k and s k are the material parameters, which are usually fitted to the experimental stressstrain curve. The initial slip resistance is assumed the same for all slip systems and is denoted s cr;0 . The loading-unloading conditions of the regularized porous crystal plasticity model are defined in Kuhn-Tucker form as and used to determine the plastic parameter _ k.

Temporal integration algorithm
In the following, a stress-update algorithm is devised for explicit finite element simulations (Hallquist 2006), assuming that the time increments are small. The stress state and all state variables are updated by the cutting plane algorithm (CPA). It is assumed that all quantities at time t n are known, e.g. F n ,r n , C n , s a cr;n , f n and R n , and in addition the deformation gradient F nþ1 at time t nþ1 is known. The rotation tensor R 0 is initialized at the start of the simulation using the initial Euler angles of the crystal and the slip resistances are all given the same initial value, i.e., s a cr;0 ¼ s cr;0 . The initial value of the void volume fraction is f 0 .
The velocity gradient L at time t nþ1=2 ¼ ð1=2Þ t n þ t nþ1 ð Þ , i.e., at the half-step, is estimated as: where Dt nþ1 ¼ t nþ1 À t n is the time increment, and thus the rate-of-deformation and spin tensors are given by: The rate-of-deformation tensor is transformed to the co-rotational coordinate system according to: It is useful to introduce the incremental strain and plastic strain tensors: and, analogously, the incremental plastic rotation tensor: First, the incremental plastic strain and rotation tensors Dê p nþ1 and Dx p nþ1 , the slip increments Dc a nþ1 ¼ _ c a nþ1=2 Dt nþ1 ; and the iterative change of the plastic multiplier dk are initialized to zero. The trial stress tensor is defined by: With the trial stress tensor calculated, the iteration scheme of the CPA is initialized. The value of the yield function with the current stress is evaluated. To this end, the stress tensorr i is first used to calculate the resolved shear stress s a i according to Eq. (11), where i is an iteration counter. Then, after the von Mises stress r vM;i and hydrostatic stress r H;i are calculated, Eq. (12) modified according to Eq. (13) is solved to obtain the effective resolved shear stress s a eff;i . The latter quantity is then used together with the slip resistances s a cr;i to calculate the value of the yield function U i . In the first iteration,r i ¼r trial i ¼ 0 ð Þ. If the value of U i is less or equal to a tolerance parameter (tol ¼ 1:0 Â 10 À10 ), then the stress is lying within the yield surface and the iterations are stopped with the current stress value. If U i is larger than the tolerance parameter, then the stress is outside the yield surface and the return map is initiated to re-establish consistency.
To this end, the yield function is linearized about the current state of the material and the result is: where the partial derivatives of U and the iterative changes of the independent variables are given in Appendices A and B. Based on this linearization, the iterative change of the plastic multiplier is calculated as: where the minor symmetry of the elasticity tensor was used, and the auxiliary variable H i is introduced as With dk iþ1 given, the incremental plastic strain and rotation tensors and the plastic slip increments are updated according to: Dx Dc The accumulated total slip is updated by: Using Eq. (21) the updated slip resistances s a cr;iþ1 are found as: The void volume fraction is updated by: and finally, the stress tensor is updated by the plastic corrector: At this point, the iteration number i is incremented and the iteration is repeated.
When convergence is reached and the magnitude of U iþ1 is below the tolerance value, the iterations are stopped. The final values of the stress tensorr nþ1 , the plastic rate-of-deformation tensorD p nþ1=2 , the plastic spin tensorŴ p nþ1=2 , the slip resistance s a cr;nþ1 , the accumulated slip C nþ1 and the slip rates _ c a nþ1=2 are obtained from the values of the last iteration. The plastic spin tensor is then rotated to the global coordinate system: To update the rotation tensor, first the incremental elastic rotation Dx e nþ1 is found as: The rotation tensor is updated using the second order update: Finally, the stress tensor is rotated back to the global coordinate system: The crystal orientation evolution may then be analysed by extracting the updated Euler angles from the rotation tensor.
In finite element simulations, the element is eroded when the void volume fraction reaches a critical value, f max . The magnitude of the effective resolved shear stress s a eff generally increases with increasing void volume fraction. From Eq. (12) it may be seen that for very high values of s a eff all but the last two terms approach zero, and the equation degenerates to 1 À q 2 1 f 2 ¼ 0 ð48Þ In the simulations, this is reflected by the asymptotical growth of s a eff when f approaches the value of q À1 1 . The elements are thus deleted when the void volume fraction is close to q À1 1 and the numerical instabilities caused by the asymptotic growth of s a eff are detected, i.e., the high values of s a eff =s a cr start to produce NaN type values in the yield function calculations. The value of the critical void volume fraction in the present model should not be confused with the critical porosity at coalescence, f c , often used in fracture studies (see e.g. Frodal et al. (2020)). In this study, the accelerated void growth, which is commonly included into Gurson-type models to account for coalescence, was not implemented to keep the number of model parameters as low as possible.

Numerical study
The implemented porous crystal plasticity model was tested for two cases. In the first case, a 2D model of a polycrystal is subjected to plane-strain tension, whereas in the second case, a uniaxial tension test of a 3D polycrystal is simulated.
The square model of a 2D polycrystal consisted of 384,400 plane-strain elements with reduced integration and Flanagan-Belytschko stiffness-based hourglass control (Flanagan and Belytschko 1981). The explicit solver of the nonlinear finite element code LS-DYNA (Hallquist 2006) was used in the calculations. Mass-scaling was applied to reduce the computational time and the kinetic energy was controlled at every step to ensure that it was very small compared to the total energy, ensuring that the simulation remained quasi-static. The solution converged to the same values for both the stress and strain fields and the global force for various time steps tested.
The Euler angles assigned to the 384,400 elements were taken directly from the calculated Euler angles of 384,400 grid points of an electron backscatter diffraction (EBSD) scan of an AA7075-T651 alloy specimen presented in Fourmeau et al. (2015). The results of the scan are presented in Fig. 1, and the distribution of one of the Euler angles (u 1 ) in the finite element model is presented in Fig. 2. The EBSD data is inherently noisy, which may be seen in Fig. 2, whereas the plot in Fig. 1 is smoothed. The polycrystal finite element model in Fig. 2 is surrounded by a layer of elements governed by von Mises plasticity with isotropic hardening defined by the average stress-strain curve of the material taken from Fourmeau et al. (2015). The hardening parameters of the crystals are obtained by a fitting procedure using the LS-OPT software, as described e.g. in Khadyko et al. (2017). The material parameters used in the simulation are summed up in Table 1. The left edge of the model is fixed, and the velocity of the right edge is ramped up to a constant value to simulate plane-strain tension. The isotropic plasticity elements provide a more natural and softer boundary for the polycrystal than fixed edges or periodicity as boundary conditions. The porous plasticity model parameters are partly taken from Han et al. (2013) and summarized in Table 2. The critical void volume fraction was taken as f max ¼ 1=q 1 % 0:67. Element erosion is used to describe crack propagation and the element is eroded when f ¼ f max in the single integration point of the finite element or when the aspect ratio (i.e., the ratio between the longest and the shortest diagonal) of an element became greater than 10. The latter criterion was used to remove elements that were heavily deformed without developing sufficiently high porosity for fracture to occur.
Plots of the von Mises equivalent plastic strain (defined as the time-integrated von Mises norm of the plastic strain rate tensorD p ) are presented in Fig. 3.
Initially the polycrystal deforms plastically with a strong tendency to form shear bands, which may start as smaller local bands and then coalesce into larger bands shearing through multiple grains. The first elements reach the critical void volume fraction f max in the point of intersecting shear bands. After that, the crack propagates along one of the shear bands.
The simulation was run utilizing eight cores of a computer cluster node (Intel Xeon X5680). The fracture initiation was reached relatively fast, in approximately 20 h. After that point though, the simulation slowed down significantly and took approximately 4 days until the specimen was separated in two halves. At the point of fracture initiation, all elements were eroded due to void growth, i.e.,  f ¼ f max , but as the simulation progressed the proportion of elements eroded due to high aspect ratio increased and a majority of the elements eroded at the later stages of fracture were removed due to a high aspect ratio. Overall, almost 40% of all elements was eroded by the aspect ratio criterion. This criterion is further discussed in Sect. 5. The second test case is a 3D model of a polycrystalline specimen with rectangular cross-section subjected to uniaxial tension along ED. The cross-section thickness to width ratio is 1:3, while the length of the polycrystal along the tensile axis is equal to two times the thickness. The polycrystal structure was generated as a Voronoi tessellation with 3000 equiaxed grains. The model consists of 750,000 8-node brick elements with reduced integration and Flanagan-Belytschko stiffness-based hourglass control. Each grain is represented on average by 250 elements and there are on average 8 grains through the thickness direction of the model. The grains are approximately equiaxed due to isotropic seeding and propagation of the Voronoi tessellation. The chosen set-up provides relatively coarse realization of the grains and jagged ''staircase''-like grain boundaries, but the goal was here to simulate the tensile behaviour of a more realistic polycrystalline sample with a multitude of grains in the thickness direction and realistic boundary conditions instead of a smaller scale, high resolution model of a partial microstructure. The description of the intragranular fracture propagation could be improved with a high-resolution model, but the effects of localization on the fracture process would require an unfeasibly large model.
The model is presented in Fig. 4. The Euler angle data for the grains is taken from Khadyko et al. (2017). The texture is a Cube texture with a minor Goss component, typical for recrystallized aluminium Fig. 3 Contour plots of the von Mises equivalent plastic strain at a 14.1%, b 17.5%, c 21.6% global logarithmic strain and d after fracture obtained in the finite element simulation of planestrain tension Fig. 4 Finite element model of a uniaxial tension test of a 3D polycrystal alloys. The material parameters governing elasticity, yielding and work-hardening are the same as in the first case, see Table 1, and the same holds for the parameters governing porous plasticity, see Table 2. Considering that we were mostly interested in testing the regularized porous crystal plasticity model and observing trends, rather than trying to quantitatively model the experimental data, we decided to use the same material parameters for the whole study. The said parameters describe an alloy with high yield strength and low work hardening, providing an early strong localization, which is more illustrative. The left and right sides of the model are again connected to parts with von Mises plasticity, which are in turn fixed on the left edge and displaced with a constant velocity on the right edge. This provides more realistic boundary conditions and helps to initiate necking in the middle of the model. Element erosion is applied to model crack propagation when f ¼ f max in the single integration point of the hexahedral elements.
The results of the simulation are presented in Fig. 5. In this case, the distinct sharp shear bands do not form, because of the less constrained plastic flow, and instead diffuse necking in the thickness direction is observed. Fracture initiates at one location on the right side of the specimen and progresses towards the left side, capturing in a realistic way the crack propagation process in polycrystalline materials under tension. The fracture surface in Fig. 6a may be compared to the fracture surface of an AA6063 aluminium alloy with the same texture (but different strength and work hardening characteristics) subjected to uniaxial tension using flat specimens with the same width to thickness ratio, obtained in Khadyko et al. (2019) and presented in Fig. 6b. The overall proportions of the fracture surface for the specimen lying in the ED direction were reproduced quite well. The simulation was run on the same eight-core node of a cluster and took about 6 days. In the 3D simulation, there were only some few elements that were eroded due to an extreme aspect ratio.
The global nominal stress-strain curves for the simulations are presented in Fig. 7. The reduction of the nominal stress (force) in these curves reflect both the strain localization (necking and shear banding), texture evolution and damage softening. The fracture surface and force reduction could be obtained also by the phenomenological plasticity models but using the crystal plasticity model allows for including texture evolution effects and microstructure heterogeneity in a natural way. The nominal stress for the plane-strain tension simulation does not fall to zero due to the elements along the boundary governed by von Mises plasticity that were not eroded. The differences in stress level and global failure strain are mostly caused by the higher constraint level in plane-strain tension than in uniaxial tension and the subsequent difference in localization type, but may also be influenced by the Fig. 5 Contour plots of the von Mises equivalent plastic strain for the 3D finite element model of the uniaxial tensile test of the 3D polycrystal a before fracture at 41.2% global logarithmic strain, b during the fracture process, and c after full fracture difference in the crystallographic texture applied in the two cases.

Discussion
The longstanding problem of rate-independent crystal plasticity is the Taylor ambiguity, or a situation where several sets of active slip systems are equally valid decompositions of the plastic strain rate. This problem is most often circumvented by using a power law type of viscoplasticity, which in geometric terms replaces a sharp vertex of the yield surface polytope with a smooth curved ''vertex''. By choosing a very low value of the rate sensitivity parameter, rate-independent behaviour can be closely approximated (Mánik and Holmedal 2014). In Paux et al. (2015), another step was made and the power law was used to combine all slip systems in a single regularized yield function with a single plastic multiplier. The variation of yield function regularization for crystals proposed in Zamiri and Pourboghrat (2010) uses the KS-function (Kreisselmeier and Steinhauser 1980), instead of a power law. In all these models the sharp vertex is replaced by a smooth one, which means that if the stress is lying exactly in the vertex of the yield locus, producing the same resolved shear stress on adjacent slip systems, it will activate all of these slip systems simultaneously, unlike the rate independent Schmid's law, where a choice has to be made. Thus, these models avoid the Taylor ambiguity altogether. According to the results in Zamiri and Pourboghrat (2010), the KS-function regularization provides both less computational time and better stability even for large values of the closeness parameter q, than the equivalent viscoplastic rate insensitive formulation.
The authors tried using the more usual hyperelastic viscoplastic rate insensitive CP model with the Gurson type damage but encountered excessive substepping and non-convergence for many cases of fracture simulations of polycrystals. This was the motivation for creating the present model. In addition, in the test simulations the hypoelastic model was at least 10% faster than the equivalent hyperelastic model implementation that the authors used previously.
The porous plasticity law from Han et al. (2013) has, as already mentioned, no explicit analytical solution for the effective resolved shear stress. To avoid the need for an iterative solution method, the hyperbolic cosine term was replaced here with polynomial terms of the Taylor series. An advantage of the polynomial solution is that it requires no initial condition or starting point for iterations, which may cause convergence problems that must be identified and mended for the iterative solution method to converge. The hyperbolic cosine function is approximated very well by the polynomial, producing very accurate solutions of Eq. (12), even with just 4 terms. The approximation is illustrated by plots of the hyperbolic cosine and the truncated Taylor series, presented in Fig. 8.
As may be seen from the figure, the error is increasing at higher values of the argument q 2 ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ð3=20Þ p ðr H =s a eff Þ of the hyperbolic cosine function in Eq. (12). In the extreme case when the von Mises stress and resolved shear stress are assumed to be negligibly small, compared to the hydrostatic stress in Eq. (12), the hyperbolic cosine term is limited by ð1 þ q 2 1 f 2 Þ=2q 1 f , which may become quite large for small values of void volume fraction f . On the other hand, at small values of the void volume fraction, the equation should degenerate to s a eff ¼ s a . This limits the values of the hyperbolic cosine term and correspondingly the error. The effective resolved shear stress found from the polynomial equation for some combinations of void volume fraction, von Mises stress and hydrostatic stress was compared to the iterative solution of the original equation and it was found that for realistic combinations, the solutions differed by less than 2%. The error increased only beyond that for some special cases with extremely high hydrostatic stresses and small void volume fractions. To further investigate this error, some smaller test simulations were performed with both polynomial and iterative solution of Eq. (12). The stress, plastic strain and porosity fields produced by the two methods were practically identical. However, the simulations produced only a limited set of stress-strain histories for the elements and situations might occur for which the errors are significant. Several variations of the equation for the effective resolved shear stress exist, derived by different methods, as described in Sect. 1. Also, Eq. (12) contains free parameters fitted to unit cell simulation results, and the error introduced by the polynomial expansion is thus not considered significant.
Some numerical aspects of the fracture simulations with crystal plasticity require further investigation. The mesh sensitivity is a known issue in finite element simulations even for phenomenological porous plasticity simulations and is not considered in this work. In some elements in the 2D simulation, the crystallographic orientation and the stress-strain history produced small hydrostatic stresses and consequently the Fig. 8 Approximation of the hyperbolic cosine function with the truncated Taylor series f n x ð Þ, where n is the highest power of the polynomial term evolution of the void volume fraction was slow. These elements continued to deform to extremely large strains and were artificially eroded using the arbitrarily chosen critical aspect ratio criterion before they could develop the critical void volume fraction. The porosity was evolving almost exclusively in the shear band, and practically all the elements deleted by the aspect ratio criterion still developed significant void volume fraction above 10%. As a possible solution, coalescence and accelerated void growth could have been adopted in the simulation, and then due to the high porosity in the heavily deformed elements, it is likely that the elements would be eroded due to ductile void growth instead of a poor aspect ratio. Another possibility is to introduce remeshing to improve the aspect ratio of the deformed elements. Nonetheless, the effect of the adopted critical aspect ratio on the results of the simulation should be investigated further. Whether the observed extreme deformation of the elements is a physical aspect of the fracture process or a numerical artefact remains an open question as well. However, it was found that the number of elements eroded due to the aspect ratio criterion was just a few in the 3D simulation of uniaxial tension and thus much less than in the 2D simulation of plane-strain tension.
The yield function developed in Han et al. (2013) has certain limitations. It is derived for spherical voids and does not account for the void shape evolution explicitly. Nevertheless, it provides a good approximation of the plastic behaviour of unit cell models for spherical voids in single crystals. Its implementation allows studying complex processes associated with fracture: various strain localization modes (necking, shear banding) interacting with damage, softening and fracture in an anisotropic crystal grain. The efficiency of the formulation allows for creating large and detailed high-resolution polycrystal models where these phenomena occur, as demonstrated in Sect. 4. While the model may be improved in different ways (most notable a void coalescence criterion can be added), even in its present form it can provide some interesting results in anisotropic fracture studies.

Concluding remarks
A numerical implementation of a porous crystal plasticity model in the explicit finite element method is proposed. The implementation combines a regularized rate-independent crystal plasticity formulation, based on the KS-function, with a constitutive equation for the effective resolved shear stress of the porous single crystal, which is modified here to allow for an analytical rather than a numerical iterative solution. The cutting plane algorithm is applied for the temporal integration of the rate constitutive equations and the regularized porous crystal plasticity model is implemented as a user-material model for the explicit solver of a commercial finite element code. The material model is tested for two cases of 2D and 3D polycrystals in tension, and promising results in qualitative agreement with experimental observations were obtained. The test simulations further showed that the cutting plane algorithm developed for the porous crystal plasticity model converges for all stages of deformation and fracture. In addition, reasonable computational times are obtained even for relatively large polycrystal models.