Finite element modeling of periodic polycrystalline aggregates with intergranular cracks

International Journal of Solids and Structures 90 (2016) 60–68

Contents lists available at ScienceDirect

International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr

Finite element modeling of periodic polycrystalline aggregates with intergranular cracks Nicolas Kowalski a, Laurent Delannay a,∗, Peng Yan b, Jean-François Remacle a a

Institute of Mechanics, Materials and Civil Engineering (iMMC), Université catholique de Louvain, iMMC, 4 Av. Georges Lemaître, B-1348 Louvain-la-Neuve, Belgium b Solid Mechanics Research Center, Beihang University (BUAA), Beijing, 100191, China

a r t i c l e

i n f o

Article history: Received 27 November 2015 Revised 16 March 2016 Available online 22 April 2016 Keywords: Anisotropy Internal stresses Cracked solid Texture

a b s t r a c t The present study addresses the prediction of the mechanical response of model polycrystalline aggregates in which the anisotropy of individual crystals induces high internal stresses as well as microcracks. A novel procedure is proposed in order to automatically generate finite element (FE) meshes that conform to the polycrystalline microstructure. The meshes obtained are periodic in 3-D and they contain cohesive zone elements along all grain boundaries. The FE model is applied to two materials with significant elastic anisotropy at the crystal level: stainless steel and graphite. A self-consistent equivalent inclusion solution is shown to produce a valid mean-field estimate of the effective stiffness of the polycrystal. This reference solution is used in order to optimize the assignment of lattice orientations to individual grains prior to the FE simulations. It is demonstrated that relying on such model polycrystals increases the accuracy of FE predictions at an affordable computational cost. The influence of internal cracks on the macroscopic response depends on the anisotropy of individual crystals. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction At the microscopic scale, polycrystalline materials deform heterogeneously. This is due to the presence of harder constituents (reinforcing phase or precipitates) and also simply to the anisotropy of individual crystallites or “grains”. Significant stress concentrations occur at the junction of adjacent grains so that grain boundaries and triple junctions are often regarded as preferential crack nucleation sites (Luther and Könke, 2009; McDowell and Dunne, 2010; Onck and van der Giessen, 1998; Sauzay and Jourdan, 2006; Simonovski and Cizelj, 2015; Yan et al., 2016). Computational modeling of the latter phenomena requires three ingredients: (i) a valid representation of the mechanical responses of microscopic constituents, (ii) the definition of a “digital microstructure” that is statistically representative of the real polycrystalline aggregate, and (iii) a robust and efficient numerical solver. It is a challenge to generate model polycrystals which warrant valid simulations of the material response at an affordable computational cost. Microstructures differ greatly depending on whether the material was produced by solidification, sintering, or a thermomechanical treatment involving plastic deformation and annealing. Models must hence account for various characteristics, includ∗

Corresponding author. Tel.: +3210472356; fax: +3210472180. E-mail address: [email protected] (L. Delannay).

http://dx.doi.org/10.1016/j.ijsolstr.2016.04.010 0020-7683/© 2016 Elsevier Ltd. All rights reserved.

ing the crystallographic texture (Bunge, 1982), the statistical distribution of grain shape and size and the topology, e.g. number of first neighbors (Groeber et al., 2008; Rowenhorst et al., 2010; Saylor et al., 2004). The required level of accuracy of the microstructure representation depends on the aim of the simulation. On the one hand, the overall elastic stiffness and yield stress are often properly apprehended based on simplified representations of the grain shapes, including cubes (Raabe et al., 2002; Sarma and Dawson, 1996), dodecahedra (Mika and Dawson, 1998) or truncated octahedra (Delannay et al., 2006; Ritz and Dawson, 2009; Zhao et al., 2007). On the other hand, localized phenomena, such as internal stresses, plastic stored energy, void nucleation or cracking, require more accurate estimates of the heterogeneity of microscopic mechanical fields and, hence, richer representations of the microstructure. Whereas the most realistic grain shapes and size distributions are obtained after simulated annealing processes (Cruz-Fabiano et al., 2014; Sieradzki and Madej, 2013), it is common to rely on a definition of the grains based on Voronoï (Diard et al., 2005; Dobrzynski et al., 2011; Fritzen et al., 2009; Gonzalez et al., 2014; Quey et al., 2011; Resk et al., 2009; Sweeney et al., 2012) or Laguerre–Voronoï (Hitti et al., 2012; Lavergne et al., 2013) cells. Besides, the model polycrystal must contain a sufficient number of grains having lattice orientations which sample the crystallographic texture in each phase. Indeed, since digital mi-

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

crostructures have finite size, simulations performed on two such microstructures with the same statistical properties do not lead to exactly the same prediction. Moreover, if the digital microstructure comprises an insufficient number of microscopic constituents, the outcome of the simulation is systematically biased relative to the average properties of a larger polycrystal. The latter is due to the effect of boundary conditions on the predicted mechanical fields (Kanit et al., 2003; Niezgoda et al., 2010; Ranganathan and OstojaStarzewski, 2008). A model microstructure is said to form a “representative volume element” (RVE) if it allows reproducing the “effective” asymptotic response of an infinitely large polycrystal with the desired accuracy. In order to apply a numerical solver, the volume of the model polycrystalline aggregate must be discretized and this can be done in different ways. The procedure proposed here generates a finite element (FE) mesh which conforms to all grain boundaries and which is periodic in 3-D. Such “conforming” meshes (also adopted e.g. in Fritzen et al., 2009; Gonzalez et al., 2014; Quey et al., 2011) allow direct modeling of intergranular cracks using cohesive zone elements laying precisely along the grain boundaries (Luther and Könke, 2009; Simonovski and Cizelj, 2015; Yan et al., 2016). Polycrystalline microstructures may also be modeled using nonconforming meshes. One option is to use the level set approach in order to define the grains and to immerse the microstructure in an unstructured finite element mesh (Dobrzynski et al., 2011; Resk et al., 2009). In the simplest case, the RVE is discretized in the form of an array of voxels, as is required when relying on fast Fourier transforms (Brenner et al., 2009; Eisenlohr et al., 2013; Lavergne et al., 2013). The paper outline is as follows. The proposed mesh generation technique is presented in Section 2. Then, Section 3 describes an original method aiming at a reduction of the RVE size (and hence the computational cost) based on an efficient assignment of crystal lattice orientations to individual grains. In Section 4, the FE model is applied to polycrystals in which strain heterogeneity arises due to elastic anisotropy of the grains. Unlike previous similar studies (Böhlke et al., 2010; Kamaya et al., 2007; Kanit et al., 2003), the present investigations considers the extreme anisotropy of graphite crystals (Yan et al., 2016) as well as the influence of intergranular cracks on the polycrystalline response. 2. Generation of a FE mesh conforming to a periodic polycrystalline microstructure The proposed finite element modeling applies to a variety of polycrystals but it is tested here only in the case of equiaxed grains represented by Voronoï cells. This means that we consider a number of seeds distributed inside a cube of size L × L × L and we then define the ith grain as the region that is closer to the ith seed than to any other seed. An advantages of Voronoï cells is that they are convex by construction. Grains considered here are delimited by planar boundaries and the intersection of three adjacent cells is a straight line segment. The microstructure is made periodic by replicating each seed 26 times: each replicated seed is translated at distances +L or −L along one edge of the cube. The microstructure generator has been implemented in Gmsh (Geuzaine and Remacle, 2009), using Voro3D (Dupuis et al., 2005) in order to compute the periodic Voronoï diagram. The positioning of the seeds must be controlled because it determines the grain size distribution and also because an inadequate positioning induces very small line segments at the boundaries of some of the cells. When the length of the latter segments is small compared to the average finite element size, the mesh tends to contain badly shaped finite elements and this, in turn, affects the predictions. Some arbitrary small Voronoï edges may exist when five seeds are close to being co-spherical, which turns

61

out to be frequent when more than ∼ 50 grains are randomly positioned. Such small edges can be avoided using a regularization procedure (Nygards and Gudmundson, 2002; Quey et al., 2011). In the present study, the seeds are initially  positioned randomly but with a minimum distance equal to L 3 5/(π ngr ) between them. Then, the microstructure is regularized by iteratively repositioning the seeds until the length of the smallest Voronoï edge reaches a threshold set equal to the average size of finite elements. Another solution would have been to remove small edges during a post-processing operation on the polycrystal geometry. However, Voronoï cells would then have non planar faces and the periodic structure of the mesh might be broken. Starting from a periodic Voronoï diagram, a set of topological entities that correspond to the points, edges, faces and regions of the Voronoï diagram is build. The topological adjacencies between those entities are subsequently computed to form a true boundary representation of the microstructure. The meshing process is then applied to the edges of the model to form a 1-D mesh. Onedimensional meshes are then used as input for generating the 2-D meshes of the faces of the model. Periodicity is managed using a master-slave approach: a master model entity is chosen for every entity that has a periodic counterpart. Master entities are meshes in a first step. Then, slave entities use translated copies of master meshes. The meshes generated here are uniform in size. This is clearly not a limitation and graded meshes can be generated using the same approach. Fig. 1 shows some meshes obtained for microstructures containing respectively 8, 27, 64 and 125 grains. Inside these meshes, which are generated in less than 3 s on a personal computer, each grain is discretized into about 1200 integration points. Cohesive zones are used to model inter-granular cracks (Segurado and LLorca, 2004). They are introduced along all grain interfaces. Consider a triangle t(v1 , v2 , v3 ) located at the interface between regions r1 and r2 . Let us call T1 (v1 , v2 , v3 , v4 ) and T2 (v3 , v2 , v1 , v5 ) the two tetrahedra that are adjacent to t and that belong respectively to r1 and r2 . A new triangle t  (v1 , v2 , v3 ) is created with vertices v1 , v2 , v3 located at the same positions as vertices v1 , v2 , v3 . We then modify the vertices of T2 (T2 (v3 , v2 , v1 , v5 )) and create a prism P (v1 , v2 , v3 , v1 , v2 , v3 ) of zero volume, forming a cohesive zone. This simple procedure introduces topological holes in the domain. Indeed, each Voronoï vertex is duplicated as many times as the number of its adjacent faces but cohesive relations only occur between copies that share a face. The impact of the (zero volume) holes is negligible in practice. Introducing cohesive zones at periodic interfaces requires special care. Here, we ensure that the cohesive elements that are related to periodic interfaces have the same periodicity relationships as the rest of the mesh. For each pair of boundary faces of the tessellation linked by a periodicity relationship, we select one that will be duplicated to form a cohesive zone, while the other will not. Within each pair of periodic faces, we duplicate the one for which the barycenter is closest to the origin. This ensures that the ensemble of duplicated faces forms a continuous surface and it simplifies the formulation of periodicity conditions. 3. Definition of the lattice orientations of the grain sampling using a reference mean-field solution Once the model microstructure is defined, it is necessary to assign a different lattice orientation to every grain constituting the polycrystalline aggregate. At this stage, the sampling of grain sizes is prescribed and it is non-uniform unless all grains have the same shape, i.e. unless Voronoï seeds are positioned on a regular grid. We further assume that the crystallographic texture is known and that the lattice disorientation of adjacent grains is arbitrary. Hence, the method does not account for a possible clustering of the grains,

62

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

Fig. 1. Examples of finite element meshes representing polycrystals containing 8, 27, 64 and 125 grains respectively. The coloring of each grain is according to the orientation of the crystal lattice.

nor for a minimum disorientation of neighboring grains that might result from recrystallization. The aim of the suggested procedure is to identify the sampling of crystal orientations that allows the best prediction of the effective elastic stiffness of the modeled material. The “effective” stiffness C refers to a polycrystal of infinite size whereas FE models allow computing the “apparent” stiffness  C of an RVE with finite size. Calling C (g) the stiffness of grain g, we have:

 C : ε ≡ σ =



fg C (g) : ε(g)

and

g=1..ng

 

C ≡ lim  C ng →∞

ε =



f g ε (g ) =

g=1..ng

(1)

(2)

The strain concentration tensor A(g) depends on the contrast of elastic stiffness between the inclusion and the matrix and also on the Hill tensor P . The latter is a function of the ellipsoid aspect ratio and of the matrix stiffness   C . According to the self-consistent scheme, ε(g) approximates ε (g) if the far-field strain is linked to



fg A(g) : ε∞

(3)

g=1..ng

After introducing (2) and (3) into (1), the polycrystal apparent stiffness is:

 C=





f g C (g ) : A (g ) :

f g A (g )

−1

(4)

g=1..ng

If the polycrystalline aggregate has a crystallographic texture, i.e. if crystal lattices of individual grains show preferred orientations, the expected overall mechanical response is anisotropic and the sampling of lattice orientations considered in the FE model must be statistically representative of the orientation distribution function (Böhlke et al., 2010; Bunge, 1982). The lattice orientations assigned to the grains must then be selected using a texture discretization algorithm e.g. Melchior and Delannay (2006) so that simulations capture the anisotropy of the real material. In the present study, we disregard texture-induced anisotropy and rather consider polycrystals that are isotropic due to a random distribution of lattice orientations. If the average grain shape is equi-axed, the effective stiffness of the polycrystal is isotropic (at least until the onset of intergranular cracking). Using μ and κ to denote the effective shear and bulk moduli, we have:

C = 2μ I dev + 3κ I vol dev





g=1..ng

where . and .(g) represent the average values computed, respectively, on the volume of the model polycrystal (i.e. on the RVE) and on the volume of grain g (which is a fraction fg of the total volume). For now, strain heterogeneity throughout the polycrystal is computed using a self-consistent homogenization scheme that is based on Eshelby’s solution of the isolated inclusion in a uniform infinite medium e.g. Brenner et al. (2009) and Dvorak (2012). According to Eshelby’s seminal work, an ellipsoidal inclusion embedded in a uniform infinite medium undergoes a uniform strain ε(g) that evolves linearly with the far-field strain ε∞ :

  ε(g) = [I − P : C − C (g) ]−1 : ε∞ = A(g) : ε∞

the polycrystal average strain as follows:

(5)

+ I vol is the fourth-order identity tensor, I vol =

where I = I 1 3 I  I and I is the second-order identity tensor. Since grains are represented by spherical inclusions in an isotropic matrix, the Hill tensor too is isotropic (Dvorak, 2012):

P=

3 ( κ + 2μ ) 1 I dev + I vol 5μ ( 3κ + 4μ ) 3κ + 4μ

(6)

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

The effective shear and bulk moduli are a priori unknown but they may be computed in different ways. For instance, an elegant analytical solution, that is specific to cubic crystal symmetry, was developed by Böhlke et al. (2010). Another approach is adopted here, which takes advantage of the fact that, in a single-phase polycrystal where all grains have the same anisotropic elastic properties, C (g) varies due solely to the differing lattice orientations. Hence, the isotropic projections I dev :: C (g) and I vol :: C (g) are uniform throughout the polycrystalline aggregate. The same applies to the strain concentration tensor A(g) on the condition that the polycrystal response is isotropic. As a consequence, μ and κ are efficiently derived from (4):

2μ =

I dev :: C (g) : A(g)  I dev :: A(g) 

=

I dev :: (C (g) : A(g) ) I dev :: A(g)

(7)

and

3κ =

I vol :: C (g) : A(g)  I vol :: A(g) 

=

I vol :: (C (g) : A(g) ) I vol :: A(g)

(8)

As the strain concentration tensor A(g) depends on μ and κ , the latter elastic constants must be computed iteratively. Starting from A(g) = I , the iterations proceed as follows: 1. Compute μ and κ using (7) and (8) and compare them to the outcome of the previous iteration. 2. If the difference is larger than the chosen threshold, update P using (6) and A(g) using (2). Then, move on to the next iteration. It has been checked that, even in the case of extreme anisotropy (graphite crystals), the largest relative misfit in the successive estimates of the elastic constants was smaller than 10−8 after less than 30 iterations. At this stage in the development, both the effective stiffness C and the strain concentration tensors A(g) are known. The assignment of individual grain orientation may thus proceed as follows. The ng grains forming the model polycrystals are first assigned random lattice orientations. Then, the ng lattice orientations are tuned in such a way as to minimize the discrepancy between the apparent stiffness  C and the reference effective stiffness C . Both of the latter tensors are conveniently expressed in the “macroscopic” reference frame {ei } that is rotated onto the “crystal” reference frame {eˇ i } using a proper orthogonal tensor R(g) . The grain stiffness tensor is equivalently expressed as g) g) C (g) = Ci(jkl ei  e j  ek  el = Cˇi(jkl eˇ i  eˇ j  eˇ k  eˇ l

(g )

where Ci jkl

(9)

(g ) (g ) (g ) (g ) = Rim R jn R R Cˇmnop . Indeed, the stiffness compoko

lp

nents in the crystal reference frame are identical in all grains (g ) (∀g, Cˇmnop ≡ Cˇmnop ). Recalling that f(g) denotes the volume fraction of each grain, Eq. (4) may be reformulated as:

 C =

ng 

g) (g ) f (g) Ci(jmn Amnkl ei  e j  ek  el

g=1

:

ng 

g) f (g) Ai(jkl ei  e j  ek  el

−1 (10)

g=1

(g ) (g ) (g ) (g ) (g ) where Ai jkl = Rim R jn Rko Rl p Aˇmnop . Each lattice orientation is de-

fined by three Euler angles (Bunge, 1982). The ng optimal sets of Euler angles are determined by minimizing  C − C  according to an iterative steepest descent algorithm, thus by making use of the derivatives of the residual relative to the orthogonal tensors R(g) . 4. Results Given the goals set for the study, we have selected two materials with significant elastic anisotropy at the level of indi-

63

Table 1 Effective elastic properties of stainless steel and graphite isotropic polycrystals computed based on the Voigt model, the self-consistent (SC) model and the Reuss model. All values expressed in GPa. Stainless steel

μ κ E

Graphite

Voigt

SC

Reuss

89 160 225

76 160 196

60 160 160

μ κ E

Voigt

SC

Reuss

220 286 525

54 89 135

10 36 28

vidual crystallites. Both materials are used in structural components of nuclear power plants. The first one is AISI304 stainless steel with Cˇiiii = 204.6 GPa, Cˇii j j = 137.7 GPa and Cˇi ji j = 126.2 GPa. The Zener coefficient (2Cˇi ji j )/(Cˇiiii − Cˇii j j ) = 3.8 indicates high anisotropy for cubic crystal symmetry. The second material, polycrystalline graphite, has hexagonal crystal symmetry. It is isotropic in the “basal plane” but extremely anisotropic along the normal eˇ 3 to the basal plane. The elastic constants of graphite are: Cˇ1111 = Cˇ2222 = 1060 GPa, Cˇ1122 = 180 GPa, Cˇ1133 = Cˇ2233 = 15 GPa, Cˇ3333 = 36.5 GPa and Cˇ1313 = Cˇ2323 = 4.5 GPa. The effective shear and bulk moduli μ and κ have been computed using (7) and (8). The corresponding Young modulus and Poisson ratio are listed in Table 1 where they are compared to the predictions of the Voigt model, assuming iso-deformation of all grains, and to the predictions of the Reuss model, assuming isostress. The latter two models constitute upper-bound and lowerbound estimates of the overall stiffness e.g. Dvorak (2012). The prediction of the bulk modulus of stainless steel is independent of the model, which is a consequence of cubic crystal symmetry. The upper and lower bounds are very distant in the case of graphite, due to the huge elastic anisotropy of individual grains. Whereas the self-consistent predictions provide asymptotic effective values, all FE simulations correspond to polycrystalline aggregates with finite size. The model polycrystals are represented in Fig. 1. They contain, respectively, 8, 27, 64 and 125 grains. The FE model is evaluated first by considering a polycrystal without cracks and then by analyzing the influence of intergranular cracks on the predictions. Simulations are performed using the abaqus finite element code. Each simulation is repeated 10 times, corresponding to 10 different samplings of lattice orientations.

4.1. FE predictions before the onset of cracking Fig. 2 shows the influence of the boundary conditions of the . In the FE simulation on the predicted apparent Young modulus E case of periodic boundary conditions (PBC), the displacements of all nodes lying on the mesh boundary are connected to the displacements of nodes positioned exactly at the opposite side of the mesh e.g. Dobrzynski et al. (2011). By doing so, grains on the periphery of the aggregate behave exactly as if they were embedded in the bulk of the polycrystal. In the case of “linear displacement boundary conditions” (LDBC), the displacement of the kth boundary node u(k) is prescribed directly as a function of its position x(k) and of the polycrystal average deformation: u(k ) = ε · x(k ) .  In Fig. 2, we observe that the predicted apparent Young modulus systematically overestimates the asymptotic effective value E . As expected (e.g. Kanit et al., 2003), both the discrepancy between  and E and the scatter of the predictions obtained for different E orientation samplings decrease when the size of the RVE increases. Simulations relying on LDBC are less accurate than PBC simulations, especially in the case of graphite. Reaching a desired accuracy requires a much larger RVE in the case of graphite as compared to the much less anisotropic stainless steel.

64

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

 with the refinement Fig. 3. Evolution of the predicted apparent Young modulus E of the mesh for a model polycrystal containing 64 graphite grains meshed with second-order tetrahedra. The apparent Young modulus is normalized by the effective Young modulus E that is computed using the self-consistent model.

 as a function of the size, i.e. Fig. 2. Evolution of the apparent Young modulus E number of grains, of the RVE. Predictions are based on either PBC or LDBC and they are obtained for graphite (top) and stainless steel (bottom). The apparent Young modulus is normalized by the effective Young modulus E which is computed using the self-consistent model.

The surprisingly “stiff” predictions for the graphite aggregate are partly due to an insufficient refinement of the finite element mesh. Fig. 3 illustrates the influence of the mesh refinement on the predicted apparent Young modulus in the case of 64 graphite crystals subjected to PBC. Whereas the simulations presented in Fig. 2 involve about 10 0 0 second-order tetrahedra per grain, Fig. 3 demonstrates that more refined meshes would produce lower apparent Young moduli. However, the computational cost of simulations performed on highly refined meshes is prohibitive, especially in the case of PBC. Simulations of the RVE with 125 grains, each containing 10 0 0 second-order tetrahedra, took several days on an eight-core computer running the linux kernel. In Section 3, an original procedure was suggested in order to improve the isotropy of the FE predictions based on an optimized assignment of lattice orientations to the grains forming the model polycrystal. The efficiency of the procedure is demonstrated in Fig. 4 where we compare the apparent Young moduli that are predicted when the macroscopic loading direction is changed (re between the spectively e1 , e2 and e3 ). Clearly, the difference E stiffest and the most compliant of the three predictions is reduced when relying on the optimized sampling of lattice orientations as compared to the random sampling of lattice orientations. As expected, the benefit of the isotropization procedure is less obvious

Fig. 4. Illustration of the anisotropy of the apparent Young modulus predicted on RVEs of different sizes in graphite (top) and stainless steel (bottom). The difference between the stiffest and the most compliant direction is reduced after applying the isotropisation algorithm. All simulations assume periodic boundary conditions.

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

Fig. 5. Distribution of the amplitude of the normal stress σ n (top) and the tangential shear stress σ t (bottom) acting on planar grain boundaries inside a graphite polycrystal. Amplitudes are normalized by the macroscopic uniaxial tensile stress σ .

when the size of the RVE increases. Moreover, the anisotropy of the macroscopic prediction is much larger in the case of graphite compared to stainless steel. Selecting a different type of boundary conditions (LDBC or PBC) influences the distribution of internal stresses. Fig. 5 presents the predicted amplitudes of the grain boundary stresses in polycrystals containing 8, 27, 64 and 125 grains, respectively. To obtain statistically meaningful results, each simulation is repeated 10 times using different samplings of lattice orientations. We distinguish the normal component σ n and the shear component σ t acting perpendicular and parallel to each grain boundary interface. Only graphite is considered here. Trends are similar in stainless steel. A sharper distribution of the stress amplitude is predicted with PBC as compared to LDBC. Moreover, whereas PBC applied to the smallest RVE (eight grains) provides the correct stress distribution, the accuracy of LDBC simulations increases progressively when the RVE size is raised. The amplitudes of internal stresses predicted with LDBC are thus artificially scattered due to the constraints on the mesh outer boundaries. 4.2. Simulation of polycrystals with intergranular cracks In this section, we analyze the influence of intergranular cracks on the stiffness of the polycrystalline aggregate and on the dis-

65

Fig. 6. Reduction of the apparent Young modulus once intergranular cracking has occurred along 10% of the total grain boundary area in polycrystals containing 8, 27, 64 or 125 grains subjected to PBC. The graphite polycrystals are presented on the top and the stainless steel polycrystals, on the bottom. The Young moduli are measured either along or transverse to the loading that caused internal cracks.

tribution of internal stresses. Cracking is considered to be brittle. It occurs along the grain boundary interfaces that experience the largest tensile stresses. During the pre-loading, the polycrystal undergoes uniaxial tension and cracks hence tend to align preferentially perpendicular to the macroscopic tensile axis. Cohesive zone elements are used to probe the normal stress and the shear stress on the grain boundary interfaces. Cohesive zone elements are undamaged and hence have constant stiffness until the normal stress reaches a threshold value. Once the threshold is attained, the cohesive zone element is removed, i.e. it is assigned zero stiffness (in both shear and tension) but the interface may still withstand compressive stresses. The initial stiffness is set equal to 50 times the effective Young modulus so that the opening of the cohesive zones is negligible prior to cracking. The pre-loading is terminated when 10% of the cohesive zone elements have been suppressed. The apparent Young modulus is predicted by simulating uniaxial tension both parallel and transverse to the pre-loading. Simulations are performed on RVEs with increasing sizes (8, 27 and 64 grains) and using PBC. Each simulation is repeated five times, using different samplings of lattice orientations, all of them optimized for improved isotropy (Section 3). According to Fig. 6, the

66

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

RVE. Only the smallest RVE (eight grains) did not produce similar predictions. 5. Discussion

Fig. 7. Distribution of the amplitude of the normal stress σ n (top) and the tangential shear stress σ t (bottom) acting on planar grain boundaries inside a graphite polycrystal containing 64 grains. Amplitudes are normalized by the macroscopic uniaxial tensile stress σ . Simulations are performed using PBC before (red) and after (green) cracking along 10% of the total grain boundary area. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

apparent Young modulus measured in the direction of the preliminary tensile test, i.e. perpendicular to the average orientation of open cracks, has reduced by about 25% in graphite and only 15% in stain steel. The scatter of the predictions is less when the size of the RVE increases but the average drop of the apparent Young modulus is rather independent of the RVE size. The reduction of the Young modulus measured transverse to the preliminary tensile test is much less. In stainless steel, the transverse Young modulus is almost insensitive to the cracks. Fig. 7 compares the distribution of internal stresses, respectively, before and after cracking. Again, cohesive zone elements are used to probe the normal stress and the shear stress acting on grain boundary interfaces in polycrystals containing 8, 27, 64 and 125 grains. Only the simulations corresponding to graphite and PBC are considered. Stress values are normalized by the macroscopic tensile stress. After cracking, the scatter of the interface normal stress is reduced more than the scatter of the shear component. It is noteworthy that, once again, a small RVE with only 27 grains predicts the same distribution of stresses as a much larger

This study has addressed several issues related to the generation of model microstructures representative of polycrystalline aggregates: choosing the number of grains to be modelled (size of the RVE), assigning lattice orientations, selecting the type of boundary conditions and generating the finite element mesh. At different stages, it was convenient to rely on a reference solution of the polycrystal effective response. It was indeed assumed that the “full field” prediction of a finite element model performed on polycrystals containing increasing numbers of grains should asymptotically converge to the “mean field” prediction of a selfconsistent model based on Eshelby’s solution of the isolated equivalent inclusion. There actually is no mathematical proof that equiaxed Voronoï cells embedded in a heterogeneous polycrystal of infinite size should on average behave as spherical inclusions in a homogeneous effective medium. However, such self-consistent approach has proven to be successful in a number of previous studies (Brenner et al., 2009). This has been confirmed here in the case of stainless steel, for which FE predictions on RVEs with increasing size rapidly come close to the self-consistent prediction. Besides, in graphite aggregates, the stiffness overestimation is due to the extreme anisotropy of individual crystals. Indeed, simulations performed on less anisotropic hexagonal materials lead to trends similar to those obtained here on stainless steel (with cubic crystal symmetry). The present study has focused on isotropic polycrystals but mean field models provide valid estimates of the effective response also in the case of anisotropic elastic polycrystals (Böhlke et al., 2010). Such models are however less accurate in the case of elastic(visco)plastic responses. Indeed, Eshelby’s solution is not easily adapted from linear elasticity to the non-linear regime e.g. Dvorak (2012). Throughout the study, periodic boundary conditions (PBC) have demonstrated that they allow getting rid of the artifacts observed when a uniform strain is applied to the mesh boundaries (LDBC). Predictions obtained with the latter boundary conditions may be highly scattered and biased relative to the reference effective response, especially in the case of the very anisotropic graphite. Unfortunately, relying on PBC instead of LDBC raises the computational cost. The bandwidth of the finite element solver is increased significantly when nodes lying at opposite sides of the mesh are connected. The slowdown of the computations was dramatic when considering meshes with ∼ 105 (second-order) tetrahedra. As a consequence, the level of refinement of the FE meshes was limited. It was chosen such that simulations performed on the largest RVEs did not last more than about 48 hr. According to Fig. 3, relying on more refined meshes would produce somewhat more compliant predictions in case of graphite. The results also demonstrate the importance of properly assigning lattice orientations to the grains constituting the aggregate. The predictions obtained using optimized samplings of orientations were more isotropic (and hence more accurate) than those relying on random orientation samplings. The suggested optimization method is a rapid pre-processing of the FE simulation input. It does require that a reference mean-field solution be available. In the description of the method used in order to improve the statistical representativity of the orientation sampling, Eqs. (5)–(8) are specific to isotropic polycrystals. However, extending the method to anisotropic polycrystals is possible. It first requires to compute the anisotropic Hill tensor P and the corresponding (anisotropic) strain concentration tensor A. A second difference is that the reference anisotropic self-consistent solution must be

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

computed using a large sampling of lattice orientations rather than relying on an isotropic projection of the tensors, as proposed here. When adopting PBC to predict the macroscopic elastic stiffness and the distribution of internal stresses, the results of simulations performed on a large RVE were closely approximated by collecting the results of about five simulations performed on a smaller RVE with different (optimized) orientation samplings. Even the influence of internal cracks was reproduced (Section 4.2). However the strategy of relying on repeated simulations on small RVEs may be less successful in the case of plastically deforming polycrystalline aggregates. When crystal plasticity simulations are performed on a unique orientation sampling, the minimum number of grains commonly used in order to predict the effective response is larger than 10 0 0 (e.g. Melchior and Delannay, 2006). Plastic yielding of a polycrystal depends on the crystallographic texture but also on the grain shapes and on the distribution of grain sizes. Properly accounting for the latter microstructure characteristics may impose the use of larger RVEs than those used in the present study. According to the simulations on the polycrystals with internal cracks, the reduction of the Young modulus is larger in graphite compared to stainless steel. This is most likely due to the greater anisotropy of graphite which causes a large heterogeneity of the stresses applied to grain boundary interfaces. Cracking occurs along the interfaces undergoing the highest tensile stresses. After cracking, the stress drops down to zero along these interfaces and the apparent Young modulus now scales with the stress acting on the other interfaces. The latter is lower in the most anisotropic material. 6. Conclusions An automatic method has been proposed in order to generate finite element meshes which conform to grain boundaries in a polycrystalline aggregate made of periodic Voronoï cells. The finite element meshes are periodic and they include cohesive elements along all grain boundaries. A self-consistent prediction of the effective elastic response of the polycrystal has been derived and it has been used in a novel procedure allowing an optimized assignment of lattice orientations to the grains forming the model polycrystal. Both the fulfillment of PBC and the optimization of the lattice orientation sampling contributed to reducing the scatter of the predictions relative to the reference mean field solution. Accurate predictions of the macroscopic stiffness and the distribution of internal stresses both before and after cracking required a larger RVE and a more refined mesh in the case of the extremely anisotropic graphite crystals as compared to stainless steel. Acknowledgments L.D. is mandated by the Belgian National Fund for Scientific Research (FSR-FNRS). Computing facilities were made available through contract FRFC-2.4502.05 with FSR-FNRS. The authors are grateful to Dr. J. Segurado (IMdea, Madrid) who provided a secondorder cohesive zone element. They also thank J.F.B. Payne and A. Tzelepi (National Nuclear Laboratory, UK) for fruitful discussions about the micromechanics of graphite. References Böhlke, T., Jöchen, K., Kraft, O., Löhe, D., Schulze, V., 2010. Elastic properties of polycrystalline microcomponents. Mech. Mater. 42 (1), 11–23. doi:10.1016/j. mechmat.20 09.08.0 07. Brenner, R., Lebensohn, R., Castelnau, O., 2009. Elastic anisotropy and yield surface estimates of polycrystals. Int. J. Solids Struct. 46 (16), 3018–3026. doi:10.1016/j. ijsolstr.20 09.04.0 01. Bunge, H.-J., 1982. In: Bunge, H.-J. (Ed.), Texture Analysis in Materials Science. Butterworth-Heinemann, pp. 42–46. doi:10.1016/B978- 0- 408- 10642- 9.50 0 08-8.

67

Cruz-Fabiano, A., Logé, R., Bernacki, M., 2014. Assessment of simplified 2d grain growth models from numerical experiments based on a level set framework. Comput. Mater. Sci. 92, 305–312. doi:10.1016/j.commatsci.2014.05.060. Delannay, L., Jacques, P.J., Kalidindi, S.R., 2006. Finite element modeling of crystal plasticity with grains shaped as truncated octahedrons. Int. J. Plast. 22 (10), 1879–1898. doi:10.1016/j.ijplas.2006.01.008. Diard, O., Leclercq, S., Rousselier, G., Cailletaud, G., 2005. Evaluation of finite element based analysis of 3d multicrystalline aggregates plasticity:Application to crystal plasticity model identification and the study of stress and strain fields near boundaries. Int. J. Plast. 21, 691–722. Dobrzynski, C., Melchior, M., Delannay, L., Remacle, J.-F., 2011. A mesh adaptation procedure for periodic domains. Int. J. Numer. Methods Eng. 86 (12), 1396–1412. doi:10.1002/nme.3106. Dupuis, F., Sadoc, J.-F., Jullien, R., Angelov, B., Mornon, J.-P., 2005. Voro3d: 3d Voronoï tessellations applied to protein structures. Bioinformatics 21 (8), 1715– 1716. doi:10.1093/bioinformatics/bth365. Dvorak, G., 2012. Micromechanics of composite materials. Solid Mech. Appl. 186, 1–460. doi:10.1007/978- 94- 007- 4101- 0_1. Eisenlohr, P., Diehl, M., Lebensohn, R., Roters, F., 2013. A spectral method solution to crystal elasto-viscoplasticity at finite strains. Int. J. Plast. 46, 37–53. doi:10.1016/ j.ijplas.2012.09.012. Fritzen, F., Böhlke, T., Schnack, E., 2009. Periodic three-dimensional mesh generation for crystalline aggregates based on Voronoï tessellations. Comput. Mech. 43 (5), 701–713. doi:10.10 07/s0 0466-0 08-0339-2. Geuzaine, C., Remacle, J.-F., 2009. Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 79 (11), 1309–1331. doi:10.1002/nme.2579. Gonzalez, D., Simonovski, I., Withers, P., da Fonseca, J.Q., 2014. Modelling the effect of elastic and plastic anisotropies on stresses at grain boundaries. Int. J. Plast. 61, 49–63. doi:10.1016/j.ijplas.2014.03.012. Groeber, M., Ghosh, S., Uchic, M.D., Dimiduk, D.M., 2008. A framework for automated analysis and simulation of 3d polycrystalline microstructures.: Part 1: Statistical characterization. Acta Mater. 56 (6), 1257–1273. doi:10.1016/j.actamat. 2007.11.041. Hitti, K., Laure, P., Coupez, T., Silva, L., Bernacki, M., 2012. Precise generation of complex statistical representative volume elements (rves) in a finite element context. Comput. Mater. Sci. 61, 224–238. doi:10.1016/j.commatsci.2012.04.011. Kamaya, M., Kawamura, Y., Kitamura, T., 2007. Three-dimensional local stress analysis on grain boundaries in polycrystalline material. Int. J. Solids Struct. 44 (10), 3267–3277. doi:10.1016/j.ijsolstr.2006.09.020. Kanit, T., Forest, S., Galliet, I., Mounoury, V., Jeulin, D., 2003. Determination of the size of the representative volume element for random composites: Statistical and numerical approach. Int. J. Solids Struct. 40, 3647–3679. Lavergne, F., Brenner, R., Sab, K., 2013. Effects of grain size distribution and stress heterogeneity on yield stress of polycrystals: A numerical approach. Comput. Mater. Sci. 77, 387–398. doi:10.1016/j.commatsci.2013.04.061. Luther, T., Könke, C., 2009. Polycrystal models for the analysis of intergranular crack growth in metallic materials. Eng. Fract. Mech. 76 (15), 2332–2343. doi:10.1016/ j.engfracmech.20 09.07.0 06. McDowell, D., Dunne, F., 2010. Microstructure-sensitive computational modeling of fatigue crack formation. Int. J. Fatigue 32 (9), 1521–1542. doi:10.1016/j.ijfatigue. 2010.01.003. Melchior, M., Delannay, L., 2006. A texture discretization technique adapted to polycrystalline aggregates with non-uniform grain size. Comput. Mater. Sci. 37 (4), 557–564. doi:10.1016/j.commatsci.20 05.12.0 02. Mika, D.P., Dawson, P.R., 1998. Effects of grain interaction on deformation in polycrystals. Mater. Sci. Eng.: A 257 (1), 62–76. doi:10.1016/S0921- 5093(98)00824- 7. Niezgoda, S.R., Turner, D.M., Fullwood, D.T., Kalidindi, S.R., 2010. Optimized structure based representative volume element sets reflecting the ensemble-averaged 2point statistics. Acta Mater. 58 (13), 4432–4445. doi:10.1016/j.actamat.2010.04. 041. Nygards, M., Gudmundson, P., 2002. Three-dimensional periodic Voronoï grain models and micromechanical fe-simulations of a two-phase steel. Comput. Mater. Sci. 24 (4), 513–519. doi:10.1016/S0927-0256(02)00156-8. Onck, P., van der Giessen, E., 1998. Micromechanics of creep fracture: simulation of intergranular crack growth. Comput. Mater. Sci. 13 (1–3), 90–102. doi:10.1016/ S0927-0256(98)0 0 049-4. Quey, R., Dawson, P., Barbe, F., 2011. Large-scale 3d random polycrystals for the finite element method: Generation, meshing and remeshing. Comput. Methods Appl. Mech. Eng. 200 (17–20), 1729–1745. doi:10.1016/j.cma.2011.01.002. Raabe, D., Zhao, Z., Mao, W., 2002. On the dependance of in-grain subdivision and deformation texture of aluminium on grain interaction. Acta Mater. 50, 4379–4394. Ranganathan, S.I., Ostoja-Starzewski, M., 2008. Scaling function, anisotropy and the size of rve in elastic random polycristals. J. Mech. Phys. Solids 56, 2773–2791. Resk, H., Delannay, L., Bernacki, M., Coupez, T., Logé, R., 2009. Adaptive mesh refinement and automatic remeshing in crystal plasticity finite element simulations. Model. Simul. Mater. Sci. Eng. 17 (7), 1–22. doi:10.1088/0965-0393/17/7/075012. Ritz, H., Dawson, P.R., 2009. Sensitivity to grain discretization of the simulated crystal stress distributions in fcc polycrystals. Model. Simul. Mater. Sci. Eng. 17, 1–21. Rowenhorst, D., Lewis, A., Spanos, G., 2010. Three dimensional analysis of grain topology and interface curvature in a β -titanium alloy. Acta Mater. 58, 5511–5519.

68

N. Kowalski et al. / International Journal of Solids and Structures 90 (2016) 60–68

Sarma, G., Dawson, P., 1996. Effects of interactions among crystals on the inhomogeneous deformations of polycrystals. Acta Mater. 44, 1937–1953. Sauzay, M., Jourdan, T., 2006. Polycrystalline microstructure, cubic elasticity, and nucleation of high-cycle fatigue cracks. Int. J. Fract. 141 (3–4), 431–446. doi:10. 1007/s10704- 006- 9005- x. Saylor, D., Fridy, J., El-Dasher, B., Jung, K.-Y., Rollett, A., 2004. Statistically representative three-dimensional microstructures based on orthogonal observation sections. Metallurg. Mater. Trans. A 35 (7), 1969–1979. doi:10.1007/ s11661- 004- 0146- 0. Segurado, J., LLorca, J., 2004. A new three-dimensional interface finite element to simulate fracture in composites. Int. J. Solids Struct. 41 (11–12), 2977–2993. doi:10.1016/j.ijsolstr.2004.01.007. Sieradzki, L., Madej, L., 2013. A perceptive comparison of the cellular automata and monte carlo techniques in application to static recrystallization modeling in polycrystalline materials. Comput. Mater. Sci. 67, 156–173. doi:10.1016/ j.commatsci.2012.08.047.

Simonovski, I., Cizelj, L., 2015. Cohesive zone modeling of intergranular cracking in polycrystalline aggregates. Nucl. Eng. Des. 283, 139–147. doi:10.1016/j. nucengdes.2014.09.041. Sweeney, C., McHugh, P., McGarry, J., Leen, S., 2012. Micromechanical methodology for fatigue in cardiovascular stents. Int. J. Fatigue 44, 202–216. doi:10.1016/j. ijfatigue.2012.04.022. Yan, P., Delannay, L., Payne, J., Tzelepi, A., 2016. Micromechanistic modelling of the polycrystalline response of graphite under temperature changes and irradiation. Carbon 96, 827–835. doi:10.1016/j.carbon.2015.10.019. Zhao, Z., Kuchnicki, S., Radovitzky, R., Cuitino, A., 2007. Influence of in-grain mesh resolution on the prediction of deformation textures in fcc polycrystals by crystal plasticity fem. Acta Mater. 55, 2361–2373.