Cylindrical void growth vs. grain fragmentation in FCC single crystals: CPFEM study for two types of loading conditions

The crystal plasticity finite element method (CPFEM) is used to investigate the coupling between the cylindrical void growth or collapse and grain refinement in face-centered cubic (FCC) single crystals. A 2D plane strain model with one void is used. The effect of the initial lattice orientation, similarities, and differences between stress- and strain-driven loading scenarios are explored. To this end, boundary conditions are enforced in two different ways. The first one is based on maintaining constant in-plane stress biaxiality via a dedicated truss element, while the second one is imposing a constant displacement biaxiality factor. Uniaxial and biaxial loading cases are studied. For the uniaxial loading case a special configuration, which enforces an equivalent pattern of plastic deformation in the pristine crystal, is selected in order to investigate the mutual interactions between the evolving void and the developed lattice rotation heterogeneity. Next, biaxial loading cases are considered for three crystal orientations, one of which is not symmetric with respect to loading directions. It is analysed how stress or strain biaxility factors and initial lattice orientation influence the void evolution in terms of its size and shape. Moreover, the consequences of variations in the resulting heterogeneity of lattice rotation are studied in the context of the grain refinement phenomenon accompanying the void evolution. Scenarios that may lead to more advanced grain fragmentation are identified.


Introduction
Nucleation, growth, and coalescence of intra-or intergranular micro-voids is a usual scenario by which ductile metallic materials fail (Benzerga and Leblond, 2010). Most often, micro-voids are nucleated as a result of decohesion or fracture process of second phase precipitates. Growth of those micro-defects takes place due to diffuse plastic deformation up to the onset of coalescence when strain localizes in the ligament connecting closely spaced voids. From the mechanics point of view, for 3D stress controlled axially symmetric loading (e.g. uniaxial tension process) the void coalescence is connected with the transition from axisymmetric to uniaxial straining mode, which leads to the plastic flow localization. After that moment the voids continue expansion, mostly towards each other up to final ligament failure or full impingement. On the other hand, the development of strain heterogeneity around a deforming and growing void leads to microstructure changes in the material around the void. It seems that the last aspect of the mechanics of porous metallic materials has not been fully explored, yet. On the other hand, this effect accompanying void evolution can have important consequences as concerns the grain refinement as observed by Beygelzimer (2005) who formulated a phenomenological model of grain fragmentation.
Beginning from the late 60s of the previous century a lot of studies have been done experimentally, theoretically, and numerically to understand the mechanics of void initiation, growth, and coalescence. Initially, numerical analyses were performed using the macroscopic isotropic nearly rate-independent elastic-plastic model of metallic material. Tvergaard (1982) and Koplik and Needleman (1988a) seem to be the first who applied respectively the 2D and 3D unit cell approach in this respect. Based on performed studies the Gurson yield criterion for porous metals, originally developed based on the analytical solution and micromechanical approach, has been modified and equipped with additional tunning parameters to become the widely used GTN (Gurson-Tvergaard-Needleman) criterion (Gurson, 1977;Tvergaard, 1982;Tvergaard and Needleman, 1984). Those initial studies revealed an important role of stress triaxiality in the void growth phenomenon. In the next studies, the influence of the third invariant (Lode parameter) and material anisotropy (Benzerga and Besson, 2001) was observed. However, as recently concluded by Pineau et al. (2016) there is not yet a universal theory of ductile failure. Moreover, as expressed by a recent review by Das (2021) there is no full agreement yet if the process is more strain or stress-driven, while these two mechanical fields are strictly related by a constitutive behaviour of the virgin medium and problem geometry (i.e. void shape and space distribution).
The void growth and coalescence contributing to the ductile failure is a process that usually takes place at the micro-level of polycrystalline metal, therefore it was soon understood that replacing the macroscopic plasticity model for a matrix material in numerical calculations with a more relevant continuum crystal plasticity framework may lead to a better understanding of the phenomenon. Following this observation, firstly, 2D analyses of unit cells with cylindrical voids were initiated under plane strain and in-plane strain-driven boundary conditions (O'Regan et al., 1997;Potirniche et al., 2006). Similar 3D analyses of spherical void growth and coalescence under strain-controlled boundary conditions were performed by Liu et al. (2007).
Those boundary conditions (i.e. strain-based) were motivated by the possibility of avoiding any instability in the calculations. The main conclusion of those studies was that the influence of crystal orientation is more significant for the loading cases with a small strain biaxility or triaxility and of secondary importance for higher strain biaxiality or triaxiality factors. Among mentioned studies, only Liu et al. (2007) provided some results related to microstructure evolution in the presence of voids. They are concerned with the texture evolution within the unit cell and the heterogeneity of deformation assessed by the misorientation angle with respect to the average orientation. It was concluded that the heterogeneity of lattice rotation is concentrated in the regions around the void.
The cylindrical void growth under plane strain while constant in-plane stress biaxiality factor in hexagonal close packed (HCP) crystal studied by Prasad et al. (2015). Yerra et al. (2010) analysed 3D cell with a spherical void maintaining constant stress triaxiality, however applied boundary conditions were not purely stress driven since at the same time equal values of two lateral macroscopic strains were imposed. The macroscopic stress direction was fully controlled in the 3D calculations by Needleman (2013, 2015) for Ni FCC single crystal, Ling et al. (2016) for FCC austenitic stainless steel, and by Selvarajou et al. (2019) for HCP Mg crystal. Different values of stress triaxiality and Lode parameter as well as selected crystal orientations with respect to loading axes were analysed. It was found that the value of the Lode parameter is more decisive concerning the void coalescence or collapse at a lower value of stress triaxiality. However, for certain anisotropic orientations, the Lode parameter can also have a significant effect on creep strain and porous evolution at higher stress triaxiality values (Srivastava and Needleman, 2015), which led to forming a polygonal void shape with rounded corners. Moreover, the initial crystal orientation dictates the location of maximum stress concentration. Based on such numerical studies analogues of the GTN criterion for single crystal were formulated by Han et al. (2013) and Paux et al. (2015), using the classical multisurface Schmid condition and the regularized Schmid law, respectively, as valid models for the bulk medium.
Let us also remark that the void growth and coalescence in the heterogeneous bulk medium described by the crystal plasticity constitutive model were also studied numerically. For example, bi-crystal unit cells were assumed in (Liu et al., 2010;Jeong et al., 2018;Dakshinamurthy et al., 2021) and polycrystal unit cells in (Lebensohn et al., 2013;Liu et al., 2021). Among mentioned contributions, selected results concerning heterogeneity of lattice rotation were provided in (Dakshinamurthy et al., 2021). Unless stated otherwise, a large strain rate-dependent CP formulation with the power law for slip was used in all papers recalled above.
The goal of the present research is two-fold. First, we would like to compare and analyse strain and stress-driven loading scenarios in the context of cylindrical void evolution in FCC single crystal under plane strain conditions. In particular, the effect of strain vs. stress in-plane biaxiality factor is elucidated. To our best knowledge, such direct comparison has not been performed in the literature yet. Second, the mutual interactions between the void evolution and development of lattice rotation heterogeneity, leading to grain refinement, as two competitive mechanisms of microstructure changes are explored. Such an interplay between two effects seems not to be sufficiently quantified in other research.
The paper is organized as follows. After this introductory section, we present the applied crystal plasticity model and its finite element implementation in Section 2. Section 3 is devoted to the description of the numerical model of a unit cell and the boundary conditions. The main body of the paper is included in Section 4 where the results of performed numerical studies are outlined and discussed. Their presentation is divided into two parts. The first concerns four uniaxial loading cases under a special crystal configuration, which enable the observation of a clear coupling of the void growth or collapse with the grain fragmentation into subgrains.
The second part of this section continues the analysis of the impact of two biaxiality factors on the void growth and overall stress-strain response for biaxial loading processes, as well as their relation to microstructure evolution. The paper is closed with conclusions.

Crystal plasticity constitutive theory
In this section, the key details of crystal plasticity implementation in FEM applied in the analyses are described. Model formulation and implementation follow Kucharski et al. (2014) and Frydrych and Kowalczyk-Gajewska (2018).
First, let us present the rate-dependent elastic-plastic model of the single crystal. In terms of kinematics description, the model follows classical contributions by Hill and Rice (1972);Asaro and Rice (1977); Asaro and Needleman (1985a). The deformation gradient F is multiplicatively decomposed into two parts: where F e and F p denote the elastic and plastic components, respectively. The evolution of the plastic part of the deformation gradient is governed by the equation: where the dot over the quantity denotes its material time derivative. The plastic velocity gradientL p is the sum of shears on slip systems: with unit vectors m r 0 and n r 0 denoting the r-th slip system direction and plane normal defined in the initial configuration. In FCC crystals, plastic deformation occurs along the {111} 110 family of slip systems, which contains 12 potentially active slip systems that are taken into account in the computations.
In the rate-dependent formulation, in order to calculate shear rates the power law (Asaro and Needleman, 1985b) is used:γ where v 0 is the material parameter,n is a rate-sensitivity parameter. The resolved shear stress τ r is the projection of the Mandel stress tensor M e on the direction and plane of slip: where S is the first Piola-Kirchhoff stress and τ is the Kirchhoff stress. The Mandel stress tensor is obtained using a hyper-elastic law: GPa   150  75 37.5 0.02 0.097 0.18 0.0 × 10 −3 1.4 1.4 10 0.001 Table 1: Elastic constants (C 11 , C 12 , C 44 ), initial critical shear stress (τ 0 ), and hardening model parameters (τ 1 , θ 0 , θ 1 , q, q 0 ), exponent in the power law (n) and reference shear rate (v 0 ) (Potirniche et al., 2006).
where C e = F T e F e is the right elastic Cauchy-Green tensor and is the Kirchhoff-type function of free energy density per unit volume in the reference configuration, L e is the anisotropic stiffness tensor of single crystal and E e = 1 2 (C e − 1) is the elastic Lagrangian strain tensor.
The evolution of the critical value of the resolved shear stress is governed by the exponential where h rs is the latent hardening parameter, equal 1 for self-hardening (r = s), q 0 for latent hardening (r = s) on coplanar systems (n r 0 ·n q 0 = 1), and q for latent hardening on non-coplanar systems (n r 0 · n q 0 = 1) and The parameters of the hardening model, elastic constants of the material, and the value of n used are shown in the table 1. The latent hardening parameter on both coplanar and non-coplanar systems is the same, but in general, it could have been taken as different.

FE implementation
The standard procedures developed for the FE implementation of finite strain elasto-plasticity in the fully Lagrangian displacement-based setting are followed (Simo and Hughes, 1998). In particular, incremental constitutive equations have been obtained by applying the implicit backward-Euler time integration scheme and the relation (2) is integrated using the exponential map, The implementation has been performed using AceGen code generator (Korelc, 2002). It combines the symbolic algebra capabilities of Wolfram Mathematica with automatic differentiation and advanced techniques of expression optimization. The package enables straightforward derivation of an algorithmic consistent tangent that leads to a quadratic convergence rate. Computations were performed using the AceFEM package. In calculations, 4-noded linear quadrilateral elements with 4 integration points are used. Additionally, the F-bar method (de Souza Neto et al., 1996) is applied in order to have a robust implementation, enabling the enforcement of nearly incompressible material behaviour in the geometrically non-linear regime. In spite of considering the 2D plane strain problem, the three-dimensional nature of the crystal plasticity model, and specifically the geometry of slip systems, are fully taken into account. At each Gauss point the material displacement gradient H = F − I is assumed for which components H i3 = H 3i = 0 (i = 1, 2, 3) and '3' denotes the direction perpendicular to the plane.

Cell model
A 2D plane strain unit cell with one cylindrical void is employed. Cartesian coordinate system (x − y) is used and the origin of the coordinate frame is placed at the node OP. The initial diameter of the void is D and the square plate has a side length of L. The ratio of D/L is used to define the void volume fraction: f = π/4(D/L) 2 . Nodes OP, XP, and YP are used to prescribe the boundary conditions. Instead of confining the sides of the unit cell to stay planar, which can over-constrain the model leading to the development of high stresses for some orientations, the periodic boundary conditions are applied. Accordingly, the displacement of corresponding nodes on opposite sides of the unit cell in the x − y plane are connected by periodic boundary conditions, namely whereH is the overall (averaged) material displacement gradient tensor of the unit cell, u 1 , u 2 represent the displacement of corresponding nodes on opposite sides and x 1 , x 2 represent the corresponding nodal vectors at the reference configuration.
Several studies were carried out by Han et al. (2013); Needleman (2013, 2015); Koplik and Needleman (1988b) on unit cells containing voids imposed with constant stress triaxiality ratio (ratio of the mean stress to the von Mises stress) boundary conditions. On the other hand, strain-controlled boundary conditions were considered by Schacht et al. (2003); Potirniche et al. (2006). In the present study for overall loading, both kinds of boundary conditions are employed in order to compare and quantify the influence of strain and stress biaxiality ratio on the void growth and coalescence. The way in which they are imposed is described in the next subsections.

Displacement controlled boundary conditions
For the displacement controlled boundary conditions, a displacement biaxiality factor β is set, which is defined as the ratio of the displacement in the x direction to the displacement in the y direction, namely β = u x (XP)/u y (YP) = const. Therefore, the following displacement boundary conditions are imposed at the reference configuration as shown in Figure 1a: • at node OP, u x = u y = 0, • at node XP, u x = βu(t), u y = 0, • at node YP, u x = 0, u y = u(t), which result in the following components of the displacement gradientH in Eq. (12) Note that all components ofH are known for this loading scenario.
For the uniaxial tension/compression case in the y-direction of the sample, considered at the beginning of the next section, the following displacement boundary conditions are imposed: • at node OP, u x = u y = 0 • at node XP, u y = 0 • at node YP, u y = u(t), which result in the following components of the displacement gradientH in Eq. (12) whereby we denoted unknown components ofH. By energy minimization, this leads to the averaged Cauchy stress for which Σ xy = Σ xx = 0, so the stress biaxiality factor η = Σ xx /Σ yy = 0.
Note that this case is not equivalent to the 3D uniaxial tension case since, in general all Σ kz , k = x, y, z are not necessarily zero for anisotropic material under plane strain conditions.

In-plane stress controlled boundary conditions
For controlling the in-plane stress biaxiality factor η, the formulation based on the proposal by Ling et al. (2016) is employed in the present study. A special spring element oriented in the direction of principal loading is employed to regulate displacement at the nodes XP and YP in order to maintain the constant stress ratio. Node OP is fixed and the displacement of node XP in the y direction is disabled to remove rigid motion as shown in the figure 1b. In-plane stress biaxiality η, which is defined as the ratio of the Cauchy's stress normal components along x direction to y direction, namely η = Σ xx /Σ yy = const is kept constant to study the void growth.
Application of the element results in the averaged Cauchy stress for the unit cell of the form, whereby we denote unknown components of Σ. Note that Σ yy (t) is also unknown, while via the truss element, displacement u y (YP) is imposed.
Different displacement and in-plane stress biaxiality ratios: β and η are considered in the present study, which are compared and discussed in the next section.

Finite element geometry and mesh
Two commercial software packages are used in this study. 2D planar model and mesh are generated using the commercial CAE software (ABAQUS version 6.13) as shown in Fig. 2.
Then the mesh data is imported to a symbolic and algebraic system, Wolfram Mathematica as specified in section 2.2 to perform finite element calculations and post-processing using the AceFEM package. The ratio of the void diameter to the side length in the x-y plane is taken as D/L = 0.2, leading to an initial void volume fraction of f = 0.0314. 2D mesh is employed  Similarly to other studies (see e.g. (Ling et al., 2016)), the overall Cauchy stress Σ = 1 JSF (J = detF) is calculated based on the volume averaged first Piola-Kirchhoff stressS found as: where V m is the bulk crystal volume and the integration is performed numerically in the reference configuration. Unknown components of the deformation gradientF = I+H are calculated based on the relation (12) using the current displacement vectors at nodes XP and YP, namelȳ

Results and Discussion
In this section, on the basis of the results of finite element simulations, the impacts of crystallographic orientation and various boundary conditions on the void growth or coalescence as well as grain refinement due to heterogeneous lattice rotation, in a 2D plane strain unit cell are examined. First, in-plane uni-axial compression and tension, simulated using the displacement controlled scenario (Eq. 13), are performed to demonstrate the void-induced heterogeneous slip activity, which then leads to spatial variation in lattice rotation. The example serves also to explore the effect of loading direction with respect to crystal axes and loading 'sign' (tension vs. compression). Moreover, the analysis preliminary verifies the model predictions with available experimental findings provided in (Gan et al., 2006). Next, various in-plane biaxial processes are considered. The displacement-controlled boundary condition will be referred to as the β loading case and the stress-controlled boundary condition will be referred to as the η loading case throughout the discussion. To see how plastic anisotropy impacts void growth in an FCC single crystal, four initial orientations of the crystalline lattice with respect to the sample axes are taken into consideration. They are collected in Table 2. In the current study, seven loading scenarios with β equal to -0.5, 0, and 1 as well as η equal to -0.5, 0, 0.8, and 1 are analysed.
The scenario η = 0.8 is selected due to its approximate equivalence to the β = 0 case. Note that the state of in-plane uni-axial tension or compression is represented by η = 0.   direction with a cylindrical void axis along the [110] (orientation O in Table 2). As discussed in detail by Gan et al. (2006), by applying the anisotropic rigid-plastic slip line theory, this configuration ensures a plane strain condition in the [001] - [110] crystal plane under the action of the compressive or tensile loading with three effective in-plane slip systems. For the pristine crystal they are results of equal activity of four systems (see Table 3), which act in opposite directions (i.e. m and −m) for tensile and compressive loading in the plane. It could be also verified that when the direction of loading is changed to [001], under a plane strain condition, the same set of slip systems will be active, again in the opposite sense. Thus, as far as plastic deformation by dislocation motion is considered, in-plane compression (cor. tension) in [001] is equivalent to in-plane tension (cor. compression) in [110].
For the sample with a cylindrical void, under the same loading conditions, the formation of regions of unequal slip activity of potentially active systems around the void is observed, which leads to the lattice rotation heterogeneity and crystal fragmentation into subgrains.
These theoretical predictions were verified by Gan et al. (2006) experimentally, for the [001] compression case, by EBSD measurements. Note that for the crystal without the void, no lattice rotation is predicted by the model, and slip activity is homogeneous, so the grain is not fragmented.
The purpose of the study is to investigate the differences between void evolution and grain fragmentation using four loading scenarios, namely tension/compression in the [110] and [001] directions, even though the same active slip systems are expected for all cases in a pristine rigid-plastic crystal (as indicated in Table 3). It is obvious that the overall stress biaxiality factor η is equal to 0 in each case.
First, we have used this example to verify the predictive capabilities of the present numerical model. As shown in Gan et al. (2006) and confirmed in our study, one of three effective inplane slip systems dominates in three different angular slip sectors which are centred at the middle of the void, as marked by dashed lines in Fig. 3b. This results in different lattice Figure 3: a) In-plane lattice rotation angle obtained using EBSD measurement (reprinted from (Gan et al., 2006) with permission from Elsevier) and b) misorientation angle map plotted found by the CPFEM method for a 5% compression strain. The dashed line represents the slip sectors at an angle of 35.3, 54.7, and 90 degrees respectively.
rotations in respective domains. In figure 3b, we present misorientation angle distribution 1 for the compression strain of 5%. Qualitatively similar subdivision is seen in experimental data quoted in figure 3 after Gan et al. (2006). There are some differences concerning the direction of rotation in the lateral domains, however, a full quantitative comparison is not possible due to the lack of the detailed experiment geometry and boundary condition data in (Gan et al., 2006).
Next, the same sample configuration is used to explore the effect of loading direction and its sign (i.e. in plane tension vs. compression) on the void evolution and grain refinement. To this end, four processes enlisted in Table 3 are studied numerically. In figure 4a, we compare the overall in-plane mean stress variation vs. magnitude of the true strain in the loading direction. It is seen that initially, the response in terms of the magnitude of the in-plane mean stress (σ mean = 1/2(Σ xx + Σ yy )) is the same for all processes and does not show visible tension-compression asymmetry. However, as the deformation proceeds, the difference starts to increase due to differences in the lattice rotation and void evolution. In each case, the stress level is smaller for the porous crystal than for the pristine one. The evolution of the normalized  decrease the void collapse is postponed to higher strain values. It should be stressed that the overall stress triaxiality value, calculated accounting for the 3D character of the stress field (note that Σ zz is not zero for all analysed cases), is approximately equal to 0.5 for both tension loadings and -0.5 for compression. Only small variations in the Lode parameter calculated for the overall stress are detected for four loading cases (its value is around 0.32-0.33).
Displacement biaxiality ratio β, calculated here as the inverse ratio of in-plane displacements in loading direction with respect to the lateral one, is seen in Fig. 4c, while the in-plane true strain biaxiality β log , calculated as the corresponding ratio of in-plane components of true strain measure (e.g. for tension/compression in [001] it is shown in Fig. 4d. Their variation with strain is compared for all four processes and pristine and voided crystals. As expected, it is observed that for a crystal without a void for all processes the evolution of β log is the same: it starts with the value of -0.5 in the elastic regime and reaches -1.0 for well-developed plastic flow, which marks incompressible deformation in that regime.
On the other hand, for voided crystals, the value of -1 is approached only for tension in [110], which is related to the stabilization of the void growth. For the remaining processes, the value does not drop below -0.95 indicating the compressibility of voided crystal.
Differences in the void growth or closing for two loading directions concern also the developed  The analysis showcased in this section illustrates that both in-plane stress biaxiality and stress triaxiality, as well as displacement or strain biaxiality, alone are inadequate in determining the growth of voids. This is particularly true when anisotropic materials are analysed. It is important to note that microstructure evolution plays a substantial role in this process.
Fragmentation of bulk crystal surrounding the void into subgrains may lead to significant impediment of the void volume changes.

Void growth and microstructure evolution in in-plane biaxial loading processes
In this subsection, to further explore factors differentiating the void growth and accompanying grain fragmentation in FCC crystals, in-plane biaxial processes are considered, for three orientations A, B, and C defined in Table 2. Orientations were selected following Potirniche et al. (2006). In order to investigate and differentiate the effect of stress and strain biaxiality seven loading scenarios with β equal to -0.5, 0, and 1 as well as η equal to -0.5, 0, 0.8, and 1 are analysed. Let us remark that orientation A, contrary to B and C, is non-symmetric with respect to the loading axes, thus shear strain component E xy (cor. shear stress component Σ xy ) may be observed for η (cor. β) loading cases even for a pristine crystal sample.

Overall response of voided crystal
Stress biaxiality ratio. The stress biaxiality ratio for seven loading scenarios is shown in Figure 7a. To start with, as it is evident, the stress biaxiality ratio for the η loading case is maintained constant for all crystal orientations during the deformation process, which verifies the validity of the finite element procedure used for imposing a constant stress biaxiality ratio.
On the other hand, in general, for displacement controlled processes (with constant β) stress biaxiality ratio η changes during the deformation process. For crystal orientations A and C, under the β = −0.5 loading case, the stress biaxiality is larger than zero; the slope initially rises, then progressively drops, and ultimately approaches the uniaxial loading case at the end of loading. Although the biaxiality ratio is marginally more than η = 0 for orientation B at the  loading. For this case, on average the value of stress biaxiality for three orientations is close to η = 0.8 that is why for comparison purposes such stress controlled scenario is also selected for analysis. Finally, for the β = 1 case, the stress biaxiality ratio is kept constant just as it does for the η = 1 case. Those graphs in conjunction with displacement biaxiality plots in Figs.
7b are important for analysing the growth of the void and the stress response. The softening stress response is evident in Fig. 7c when the stress biaxiality ratio increases and void growth is significant, resulting in coalescence.
Additionally, yellow lines are denoted as 'uniaxial' in Figs. 7, show the results obtained for the in-plane uniaxial tension process without the employment of a special spring element but using the displacement-controlled conditions with H described by Eq. 13. Calculations are performed for verification purposes and are in good agreement with the predictions obtained with the use of the spring element (marked as η = 0 in figures).
Displacement biaxiality ratio. The displacement biaxiality under various loading instances is depicted in Figure 7b. Similar to the situation of stress biaxiality, the displacement biaxiality ratio β is kept constant during the β-type process, which verifies the finite element procedure.
On the contrary, in general, for stress ratio controlled processes (with constant η), the displace-   Figure 7a shows that the stress biaxiality ratio is greater than 0 (positive) for β = -0.5, 0 and 1, and η =0.8 and 1 loading cases. As a result, in the initial deformation stage, a stiffer stress response is observed, followed by a softening response due to significant void expansion in the crystal, which cannot be further compensated by an increase of average stress in the bulk crystal. When the magnitude of peak stress for the different orientations is compared, orientation C has the largest peak stress, and orientation A has the lowest peak stress for β =1 loading case. Furthermore, the evolution of the overall mean stress in Fig. 7c correlates well with the displacement biaxiality ratio β shown in Fig. 7b. In particular, the higher β value the more stiff the initial response is and the sooner (in terms of the value of F 22 − 1) the peak stress is achieved for the given process. β and η are both equal to 0 and 1 are analysed, together with scenario η = 0.8 which approximately corresponds to β = 0 case as discussed before. It is evident that all loading cases exhibit the anisotropic response. Also, it is apparent that the response of the porous crystal differs substantially from that of the pristine crystal. With the exception of η = 0 (uni-axial loading condition), the pristine crystal response is purely elastic. In the case of η = 0 loading, the pristine crystal displays a stiffer response than the porous crystal for orientations A and C; however, for orientation B, the response is almost the same for both the pristine and porous crystal. The response of orientation C is the stiffest in each of the loading conditions. For loading scenarios, β = 0, 1, and η = 0.8, 1, orientation A initially displays the softest response, whereas orientation B exhibits the softest response by the end of the deformation process.
When the loading scenarios β = 0 and η = 0 are compared, the substantially higher stress biaxiality in the β = 0 loading case (refer to Fig. 7a) causes a more stiff response during an early deformation stage, followed by a softening due to significant void expansion in the crystal.
On the contrary, a monotonic stress increase is observed for the η = 0 loading scenario. Instead, as expected, the stress response for β = 0 case is close to η = 0.8 loading conditions. Due to the highest stress biaxiality, a similar response was observed for β = η = 1.
Normalized void volume fraction evolution. Figure 7d compares the evolution of the normalized void volume fraction for various loading conditions and the specified orientation. These evolution plots are in good agreement with displacement biaxiality ratio plots in Fig. 7a. For all orientations, the void growth rate increases as the displacement biaxiality ratio increases.
The void is collapsing for the η = −0.5 loading case, and this behaviour is the most pronounced in orientation C. Confirmation of the phenomenon can be seen in contour plots of accumulated shear for orientation C which are shown in Fig. 13. Under η = 0 loading case, the void grows very slowly as compared to higher stress biaxiality cases. If the displacement biaxiality ratio is less than -0.5, softening behaviour is not observed, since not much void growth is seen. The void growth increases at first with β = −0.5 but subsequently stabilizes for all orientations.
The evolution of the void volume fraction under the β and η = 1 loading scenarios correlates with the evolution of the displacement biaxiality ratio (Fig 7b). The curves of void evolution start to deviate from each other at the same moment when the value of β for η = 1 case drops below one.
Qualitative differences are observed in the curves shown in Fig. 7d for η-cases, which can be explained by the accompanied variation of displacement biaxility ratio. As it is seen, for the high stress biaxiality ratio: η = 1, initially for all three orientations the displacement biaxiality β is equal to 1, so the cell is under the conditions beneficial for the void expansion. Accordingly, at this stage, the void is growing in all directions (compare Fig. S.2f, S.3f, and S.4f in the supplementary material). However, as deformation proceeds the displacement biaxiality ratio is decreasing towards zero, which effectively slows down the void growth since its growth starts to be limited to one direction in the plane. Nevertheless, the void volume fraction is still growing on the cost of bulk crystal, and the achieved values are high. This causes a decrease of the overall in-plane mean stress as an increase of average stress in the bulk crystal is not able to compensate for the void expansion, Fig. 7c. On the contrary, for smaller stress biaxiality ratio: η = 0, the initial displacement biaxiality is negative, so even though the net change of void volume fraction is positive, in this scenario the void diameter is growing only in one direction while decreasing in the perpendicular one (compare Fig. S.2e, S.3e and S.4e in the supplementary material). For this case, as the deformation proceeds the displacement biaxiality ratio increases towards zero, which in this case leads to accelerated void growth because the reduction of void size in one of the directions is halted, while it is still growing in another one. For all three orientations for the considered deformation range, the increase in the void volume fraction is not yet sufficient to overcome the overall mean stress increase due to the strain hardening in bulk crystal. However, with increasing deformation one may expect softening which will be accompanied by an accelerated void growth rate. Interestingly, for the η = 0.8 case the former and latter scenarios of void growth are observed for orientations C and B, respectively (see also Fig. 9c). Orientation A exhibits here some limit case with the almost constant rate of void volume fraction. Figure 9 compares the evolution of the void volume fraction for different crystal orientations and the selected loading conditions. Five loading cases, the same as in mean stress response plots shown in Fig. 8 are illustrated. For η = 0 (in-plane uniaxial loading case), the anisotropic response is observed. Void growth in orientation C is the highest, followed by orientations A and B. But for β = 0 and 1 loading cases, due to relatively large strain biaxiality, the void growth is significant and the effect of crystal orientation diminishes. It has been verified that the latter observation is true also for other processes in which β value is kept constant. The same observation is reported in (Potirniche et al., 2006) under displacement controlled boundary conditions. The void growth for orientations C and B are nearly identical for η = 1 loading. value is tending to zero the coalescence is approached. It is seen that for the case β = 1 and symmetric orientations B and C the coalescence state is attained in two perpendicular bands along X and Y directions. Additionally, the void is loosing its spherical shape, more importantly for orientation B than C. For other cases the coalescence is approached mainly in X direction and this state is attained visibly sooner for orientations A and B than for orientation C. Figure   10a shows that, although the orientation effect is not seen in the normalized void evolution plots shown in Fig. 9 under the same value of displacement biaxiality, it manifests in the developed void shape and thus may influence the coalescence strain and, in general, the failure mode.

Local sample response
Local distribution of accumulated shear. First, in order to show the possible failure mode, contour plots of accumulated shear are presented at the end of the deformation process at strain level 0.3 for six considered loading scenarios in Figs. 11, 12 and 13, for three crystallographic orientations A, B, and C listed in Table 2, respectively. Since the strain level along the principal loading direction is the same for all the cases, one is able to observe relative variation in a shape change of the cell as a whole and the void for all six loading scenarios. Additionally, to illustrate the strain localization process, the contour plots are presented for  shape with rounded corners is observed. Similar behaviour was reported in (Srivastava and Needleman, 2015) at high stress triaxialities. Furthermore, for both orientations, substantial shear accumulation is seen around the void. The void rotates even further in asymmetric orientation A, but its shape is not perfectly circular or polygonal. The same rapid void growth is clearly observed in normalized void volume fractions plots for this β loading case and three orientations (refer Fig. 7d).
Now, let us move to the η = const loading scenarios.
Subfigures ( Finally, part (f) of figures 11-13 depicts the accumulated shear contour plots under the η = 1 loading scenario. Similarly to the β = 1 loading scenario, void growth is accelerated due to high stress and displacement biaxiality. The displacement biaxiality ratio plot (Fig. 7b) explains the slight deviations from the β = 1 loading condition. For orientations B and C, the void growth is rapid even at the strain level of 0.15, which is identical to the β = 1 loading situation. However, there is some deviation in the strain accumulation as compared to β = 1 case for asymmetric orientation A, which can be correlated with differences seen in the displacement biaxiality ratio curves in Fig. 7b for these cases. A polygonal void shape with rounded corners is seen for orientations B and C at the strain level of 0.3, which is identical to the β = 1 loading situation.
The void growth is slightly reduced since the displacement biaxiality ratio is less than 1 (refer to Fig. 7d). The primary difference is that in the β = 1 loading situation, void coalescence occurs in both loading directions, but in the η = 1 loading instance, void coalescence happens only in the transverse ligament. Displacement biaxiality plots (Fig. 7b) clearly show the origin of this disparity. In addition, the maximum shear accumulates around the void for all orientations.
Local distribution of lattice rotation. The influence of void evolution and loading conditions on new grain formation is now studied on the basis of contour plots of the lattice rotation angle.
We concentrate on, somewhat opposite, cases of orientation A and B, and present only selected results for orientation C.
The lattice rotation angle Ψ ∈ (0, π), presented in the plots, is defined as: where ∆R(t) is calculated based on the initial orientation tensor R(0) and current orientation tensor R(t), respectively, as Orientation tensor R(t) is constructed based on the current orientation of lattice direction a with the Miller indices [100] and lattice plane normal b with {001}, respectively. The change of their orientation during the deformation process is governed by the elastic part of the deformation gradient F e , so that a(t) = F e (t)a(0) and b(t) = F −T e (t)b(0). In each loading case, we observe rotation angle heterogeneity, which results in the development of a new microstructure. The presence of a void causes heterogeneity of strain, which results in heterogeneity of lattice rotation. However, we notice that the distribution of the rotation angle does not follow perfectly the distribution of accumulated shear Γ, as was already seen in Section 4.1 for uniaxial loading cases.
The lattice rotation angle plots for asymmetric orientation A are shown in Figure 14    observe lattice rotation for pristine crystal. For a voided crystal, for η = 0 and η = −0.5 loading cases, the heterogeneity of lattice rotation angle is very small, which is less than 5 • , following homogeneity of deformation seen in Fig. 12(d, e). Around the void, a few small domains with 5 − 10 • of lattice rotation are present. The formation of new grains around the void is observed for higher biaxiality loading cases where β = 0, 1 and η = 1, and the crystal domain in a larger distance from the void does not rotate significantly. Under β = −0.5 and η = 0.8 loadings, the response is somewhere in between the two scenarios discussed before.
In order to further illustrate the effects related to grain refinement, Figures 16 and 17 present histograms of the lattice rotation angle generated based on the data in Figures 14 and   15, respectively. The histogram plot on the left displays the entire unit cell, while the histogram plot on the right shows the area surrounding the void. For the purpose of the latter plot, we employed two layers of finite elements which surround the void (refer to Fig. 2). Different   Table 4: Mean (M) and standard deviation (SD) of misorientation angles (degrees) for two orientations under different loading scenarios colours of bars correspond to different loading conditions. When we compare the results for asymmetric orientation A and symmetric orientation B, we observe that orientation A has a substantially larger orientation spread than orientation B. This is because most of the elements rotate less than 10 degrees in orientation B. When we consider the area around the void for both orientations, the orientation spread widens significantly, especially for η = 1 and β = 1 loading cases.
In order to quantify more directly observed differences for each case mean value and standard deviation of rotation angle were calculated and collected in Table 4, which includes also the respective values found for orientation C. As expected, for all loading conditions the highest mean value is obtained for orientation A, which is connected with lack of symmetry for this configuration. Considering the results for the same orientation but different loading conditions, we see that the highest mean misorientation angle is found for the case η = −0.5 and 0 for orientation A, and for β = 1 for orientation B and C (refer to Table 4). The standard deviation is used to illustrate lattice rotation heterogeneity. High biaxiality factors η and β = 1 correspondingly showed the highest lattice rotation heterogeneity for orientation A. However, for other stress and strain biaxiality factors, the variation in lattice rotation is not significantly lower. The smallest heterogeneity is observed for the η = 0 case (around 4.5 degrees). As a result, in the presence of a void, orientation A is prone to grain refinement. For orientations B and C, the disparities in lattice rotation heterogeneity are larger. Again, β and η = 1 had the highest values. However, no heterogeneity is evident for η = −0.5 and η = 0 for orientation B, as shown by contour plots (Figure 15). This observation is not true for orientation C, which can be correlated with important strain heterogeneity for those two cases seen in Fig. 13. On overall, the magnitude of grain refinement appears to be more influenced by loading conditions in the case of symmetric orientations.

Summary and conclusions
In this paper, using the crystal plasticity theory combined with the finite element method, we have investigated the effects of initial crystallographic orientation, stress, and displacement controlled loading conditions on the void and microstructure evolution in a 2D plane strain unit cell. Uniaxial and biaxial loading cases have been studied.
For uniaxial loading cases a special configuration, which enforces an equivalent pattern of plastic deformation in the pristine crystal, has been selected in order to investigate the mutual interactions between the evolving void and the lattice rotation heterogeneity. It has been found that neither macroscopic in-plane stress biaxiality nor displacement/strain biaxiality, are sufficient to fully decide about the void growth, especially when anisotropic materials are considered, and that a significant role in this process is played by microstructure evolution.
Fragmentation of bulk crystal surrounding the void into subgrains may lead to significant disturbance of the void volume changes. Note that a similar observation, about the importance of the microstructure changes, was made by Prasad et al. (2015) for HCP crystal in which the appearance of domains with new twin related orientation strongly affected void growth and coalescence.
Next, biaxial loading cases have been considered for three crystal orientations, one of which is not symmetric with respect to loading directions. It has been analysed how stress or strain biaxility factors and initial lattice orientation influence the void evolution in terms of its size and shape. Overall, seven cases with three displacement controlled loading scenarios (β = {−0.5, 0, 1}) and four stress controlled loading scenarios (η = {−0.5, 0, 0.8, 1}) have been considered. The following are the key conclusions of the study: 1. It seems that the primary driving factor for void growth and coalescence is the displacement biaxiality factor β. A clearer correlation is found between variations in displacement biaxiality ratio and normalized void volume fraction evolution plots, as well as the resulting void shape and coalescence pattern.
2. Softening stress response is evident for large displacement biaxiality factors when the stress biaxiality ratio η increases. The void volume fraction increase in such cases is significant, resulting in void coalescence. The effect of crystal orientation is then diminished.
Similar findings were reported in other studies (Potirniche et al., 2006). The coalescence is observed in both directions for displacement biaxiality β =1, but only in the transverse ligament for stress biaxiality η =1. For advanced plastic deformation, particularly at high stress and displacement biaxiality η = β = 1, voids evolve into polygonal forms. Similar findings have been reported by Srivastava and Needleman (2015).
3. For stress controlled processes the starting point can be described as a biaxial straining process, which under the void growth is approaching an uniaxial straining mode. The way by which the void growth proceeds is governed by the variation of the displacement biaxiality factor β. When initially β is positive the obtained void volume fractions are larger (softening is observed earlier), while the void growth rate will be decreasing when the uniaxial straining mode is approached. On the other hand, when initially β is negative then the obtained void volume fractions are smaller (softening is observed later), while the void growth rate will be increasing when the uniaxial straining mode is approached.
4. For lower stress η and displacement β biaxiality values, an anisotropic response is observed, and the strain-stress response is dependent on crystallographic orientation. For the lowest value of stress biaxiality η = −0.5, void closure has been observed, particularly in the non-symmetric orientation A and orientation C, as well as the formation of strain localization bands.
5. In general, the heterogeneity of plastic deformation is the largest for non-symmetric orientation A. This results in lattice rotation heterogeneity and the formation of grain fragmentation in each loading case. For other orientations heterogeneity of lattice rotation is concentrated around the void, especially for higher stress and displacement biaxiality ratios (β = {0, 1} & η = {0.8, 1}). On the other hand, for small or negative values of both biaxiality factors, void evolution, and lattice rotation heterogeneity is greatly influenced by initial crystal orientation and substantially differ for the same value of stress and strain biaxiality factor, while the grain refinement encompasses a larger crystal volume.
It should be remarked that FCC crystals usually present smaller plastic anisotropy than HCP crystals, for which different types of slip systems can be activated with substantially different values of critical shear stresses. Moreover, for many HCP metals, uniaxial twinning plays an important role. In such a situation, we may expect even more significant influence of microstructure changes on void evolution and accompanying ductile failure mode. This, together with extending the analysis to 3D spherical voids, is an interesting direction for further research.