Grain-Scale Modeling Approaches for Polycrystalline Aggregates

In polycrystalline aggregates microstructure plays an important role in the evolution of stresses and strains and consequently development of damage processes such as for example evolution of microstructurally small cracks and fatigue. Random grain shapes and sizes, combined with different crystallographic orientations, inclusions, voids and other microstructural features result in locally anisotropic behavior of the microstructure with direct influence on the damage initialization and evolution (Hussain, 1997; Hussain et al., 1993; King et al., 2008a; Miller, 1987). To account for these effects grain-scale or meso-scale models of polycrystalline aggregates are being developed and are increasingly being used.


Introduction
In polycrystalline aggregates microstructure plays an important role in the evolution of stresses and strains and consequently development of damage processes such as for example evolution of microstructurally small cracks and fatigue. Random grain shapes and sizes, combined with different crystallographic orientations, inclusions, voids and other microstructural features result in locally anisotropic behavior of the microstructure with direct influence on the damage initialization and evolution (Hussain, 1997;Hussain et al., 1993;King et al., 2008a;Miller, 1987). To account for these effects grain-scale or meso-scale models of polycrystalline aggregates are being developed and are increasingly being used.
In this chapter we present some of the most often used approaches to modeling polycrystalline aggregates, starting from more simplistic approaches and up to the most state-of-the art approaches that draw on the as-measured properties of the microstructure. The models are usually based on the finite element approach and differ by a) the level to which they account for the complex geometry of polycrystalline aggregates and b) the sophistication of the used constitutive model. In some approaches two dimensional models are used with grains approximated using simple geometrical shapes like rectangles (Bennett & McDowell, 2003;Potirniche & Daniewicz, 2003) and hexagons (Sauzay, 2007;Shabir et al., 2011). More advanced approaches employ analytical geometrical models like Voronoi tessellation in 2D (Simonovski & Cizelj, 2007;Watanabe et al., 1998) and 3D (Cailletaud et al., 2003;Diard et al., 2005;Kamaya & Itakura, 2009;Simonovski & Cizelj, 2011a). In the most advanced approaches, however, grain geometry is based on experimentally obtained geometry (Lewis & Geltmacher, 2006;Qidwai et al., 2009;Simonovski & Cizelj, 2011b) using methods such as serial sectioning or X-ray diffraction contrast tomography (DCT) (Johnson et al., 2008;Ludwig et al., 2008). These approaches are often referred to as "image-based computational modeling" and can also embed in the model measured properties such as crystallographic orientations. The acquired information is of immense value for advancing our understanding of materials and for developing advanced multiscale computational models. The rather high level of available details may render extremely complex geometries, resulting in highly challenging preparation of finite element (FE) models (Simonovski & Cizelj, 2011a) and computationally extremely demanding simulations. These two constraints have so far limited the development and use

Analytical models of structures
In material science Voronoi tessellations are extensively used to model grain geometry with the purpose of calculating the properties of polycrystalline aggregates (Cailletaud et al., 2003;Kovač & Cizelj, 2005), modeling short crack initiation and propagation or modeling intergranular stress corrosion cracks (Kamaya & Itakura, 2009;Musienko & Cailletaud, 2009;Simonovski & Cizelj, 2011b). A Voronoi tessellation is a cell structure constructed from randomly positioned points, also referred to as Poisson points. For polycrystalline aggregates we can think of these points as points where the solidification starts and then uniformly extends in all directions. A solidification front then expands until it collides with another one, thus creating a grain boundary. In geometrical terms the grain boundary is obtained by introducing lines perpendicular to lines connecting neighboring Poisson points. The result is a set of convex polygons/polyhedra, see Fig. 1. Examples of 2D Voronoi tessellations can be seen in Fig. 2. Additional details on mathematical background and applications are available in (Aurenhammer, 1991;Okabe et al., 2000).
Three dimensional Voronoi tessellations can be created using the same approach. A number of mathematical and programming packages like Matlab have the option of constructing such tessellations. Fig. 3 shows examples of 3D Voronoi tessellations generated using Qhull algorithm (Qhull code for Convex Hull, Delaunay Triangulation, Voronoi Diagram, and Halfspace Intersection about a Point, n.d.) and implemented in (Petrič, 2010).        shapes of the grains are then reconstructed from the data from the 2D layers. However, the problem with this approach is that the specimen is destroyed during the measurement procedure. In last years new experimental techniques have enabled non-destructive spatial characterization of polycrystalline aggregates. Differential aperture X-ray microscopy (Larson et al., 2002), 3D X-ray diffraction microscopy (3DXRD) (Poulsen, 2004) and X-ray diffraction contrast tomography (DCT) (Johnson et al., 2008;Ludwig et al., 2008) are examples of these procedures. Through DCT for example, grain shapes and orientations can be measured and even crack initiation and growth can be monitored (Herbig et al., 2011). Since for some of the presented cases the data has been acquired using DCT, the next section gives an overview of this technique.

X-ray diffraction contrast tomography (DCT)
DCT (King et al., 2008a;Ludwig et al., 2008) is a measurement procedure jointly developed by the European Synchrotron Radiation Facility (ESRF) and University of Manchester, Materials Performance Centre, School of Materials. It combines the X-ray diffraction imaging and image reconstruction from projections to obtain the data on the grain shapes and crystallographic  orientations. A rotating polycrystalline sample is exposed to a monochromatic X-ray wave while the projection images are recorded, Fig. 4. Since the sample is rotating, each grain will pass through Bragg diffraction alignments several times. A detector system, significantly bigger than the sample, captures low index reflections. In absence of orientation and strain gradients within the grains, the diffracted beams form 2D spots that can be treated as parallel projections of the grains' volume (King et al., 2010). The shape of each grain can then be reconstructed in 3D using algebraic reconstruction techniques (Gordon et al., 1970). The resolution of the technique is in the order of 1 µm. An example of a measured grain shape of a 400 µm diameter stainless steel wire is given in Fig. 5. The complete wire is depicted in Fig. 6. One can see that the available level of detail is very high and that the obtained geometry is extremely complex. This results in highly challenging preparation of finite element (FE) models and computationally extremely demanding simulations. These two constraints have so far limited the development and use of the image-based models. However, with suitable simplifications and parallel pre-processing (Simonovski & Cizelj, 2011a), appropriate FE models can be built in a reasonable time.

Experimental data
The experimental data used in this work is of a 400 µm diameter stainless steel wire characterized in 3D by DCT (King et al., 2008a). The data has been kindly provided by the University of Manchester, Materials Performance Centre, School of Materials and comprises of 362 grains and some 1600 grain boundaries. The data provides information on the crystallographic orientation in points of a 346 by 346 by 282 grid. The experimental data can be represented as an array of 282 slices, separated in the depth (Z) direction by 1.4 µm. For each slice, crystallographic orientation has been measured on a 346 by 346 grid with 1.4 µm distance between the points on a grid (in the X and Y direction). Voxels having the same crystallographic orientation constitute a grain.

From the measured data to the surfaces
DCT characterization of a polycrystalline aggregate results in voxel-based data. Voxel-based data is also obtained in other experimental techniques like computed tomography (CT) or magnetic resonance imaging (MRI). To obtain the shapes of individual grains their surfaces need to be reconstructed from the voxel-based data.

Treatment of holes
The original DCT data contains 'holes' in the reconstructed grains, see the left-hand-side of Fig. 6. These are typical artifacts due to the limited number of projections available for each grain and the presence of erroneous contrast (Johnson et al., 2008;Ludwig et al., 2008). Since the holes are not expected in the experimental data, a simple and efficient treatment algorithm can be used to fill the holes by grain growth: • The layers with depth of one voxel are successively added to the grains in the holes vicinity.
Each layer is defined along the border between grains and holes by inspecting the one voxel deep neighborhood of the voxels within holes. Only the hole voxels with exactly one neighboring grain are added to this grain within each layer.
• Only the hole voxels with multiple neighbors remain after the first step. These are assigned to the grains dominating their immediate (one voxel deep) neighborhood.

54
Polycrystalline Materials -Theoretical and Practical Aspects www.intechopen.com The right-hand-side of Fig. 6 displays wire data-set after the treatment above.

Surface reconstruction
Geometries of individual grains need to be reconstructed from the voxel-based data. This is usually achieved through reconstructing the surfaces of individual grains. Surface reconstruction from voxels is available in a number of commercial visualization tools. The origins of these tools can mainly be traced to the field of medical visualization. The tools were later further developed for the application to material science. In this work surfaces are reconstructed as sets of triangles with Amira package (Visage Imaging GmbH, 2010). A label is assigned to each measured point, defining to which grain this point belongs. The labels are equal to the index of the crystallographic orientation. Label 1 refers to grain 1 with crystallographic orientation index 1 and so forth.
Amira's built-in SurfaceGen tool with unconstrained smoothing option is used. This tool partitions the bounding volume into 362 grains depending on the number of different labels in the 8 vertices of a given voxel. Near the triple points between the grains and near the grain boundaries vortices of a given voxel will be distributed among several grains. In these cases the voxel is subdivided into up to 6 3 sub-vortexes to give a topologically correct representation of the implicitly defined separating surfaces (Westerhoff, 2003). If two adjacent sub-vortexes are of different grains, their common face is added to the list of boundaries between the two grains. A comprehensive explanation of the procedure is given in (Stalling et al., 1998;Westerhoff, 2003) and was later implemented in Amira. Described approach automatically increases the resolution near the triple points between the grains and near the grain boundaries where vortices of a given voxel are distributed among several labels/grains. This is especially important since stress increases at these points can be expected due to different crystallographic orientations of the adjacent grains. Further details on the implemented approach can be found in (Simonovski & Cizelj, 2011a).
The density of the triangles forming the reconstructed surfaces is limited by the resolution of the experimental data. At full resolution the number of triangles is 4 758 871, resulting in FE model with 51 211 552 finite elements. The number of triangles therefore needs to be decreased. This is done using Amira's built-in surface simplification tool (Zachow et al., 2007). The simplification decreases the details as well as resolution at the triple lines between the grains, see Fig. 7 where the number of triangles has been decreased to 30 000 (case 30K), 150 000 (case 150K) and 300 000 (case 300K).

From surfaces to FE models
The complexity of the reconstructed surfaces, together with a rather large number of grains, essentially prevents the direct use of the finite element meshing capabilities of either Amira or even professional mesh generators such as for example ABAQUS/CAE (Simulia, 2010). A framework for the automatic meshing has been therefore developed (Simonovski & Cizelj, 2011a) using the Python scripting language and ABAQUS/CAE meshing tools, which are fully accessible through Python.
The framework can be applied to both analytical models of structures (e.g. 3D Voronoi tessellations) and to data obtained from experimental techniques. In both cases the surface structures are defined by the triangle-based surfaces, bounding the volume of individual grains. In the case of voxel-based data these surfaces are reconstructed with Amira. For  analytical models such as 3D Voronoi tessellations, the spatial structure is generated by the underlying analytical model (e.g. Qhull algorithm (Qhull code for Convex Hull, Delaunay Triangulation, Voronoi Diagram, and Halfspace Intersection about a Point, n.d.) implemented in (Petrič, 2010)) and surface reconstruction is not needed.
Before starting the FE meshing procedure the surface triangles aspect ratios are checked. Surface triangle aspect ratio is defined in this work as the ratio of the circumscribed circle and the inscribed circle of a triangle. Triangles with aspect ratio of more than 1000 are removed by collapsing triangle's shortest edge, removing the triangle from the structure and updating the vertices and triangles. The procedure is performed iteratively until the worst aspect ratio is above 1000. This approach improves the FE mesh quality. The triangle-based surfaces are also checked for possible errors like intersections and corrected, if necessary, by slight displacement of appropriate triangle's corner points.
The ability of exporting the reconstructed surfaces into a standard CAD format that can be read by FE pre-processors is lacking in many visualization tools. Instead, the user is encouraged to use the tool's built in meshers, which often do not match the capabilities of dedicated FE mesh engines. Furthermore, FE pre-processors do not support export formats of the visualization tools. So there is a basic difficulty of importing the reconstructed geometry into FE pre-processors. This issue is circumvented here by developing a function for exporting the reconstructed surfaces into ACIS SAT file that can be imported into practically any FE pre-processor. The surfaces of each grain are therefore saved into ACIS SAT file with all 56 Polycrystalline Materials -Theoretical and Practical Aspects www.intechopen.com Grain-Scale Modeling Approaches for Polycrystalline Aggregates 9 the surface triangles that they contain. Next, individual ACIS SAT file is imported into ABAQUS/CAE pre-processor, assigned seeds and its surfaces are meshed using triangular FE elements of a selected size. Mesh density is user controlled. The coarsest possible mesh is defined by the triangles of the reconstructed surface. Finer meshes can be obtained by using several FE elements per each triangle of the reconstructed surface. Meshing individual surfaces is independent upon each other which provides for an efficient parallelization of the process. A surface between the adjoining grains needs to be meshed only once since the shared surface is identical for both grains. To obtain a conformal mesh between adjoining grains the meshed surfaces are imprinted onto the corresponding grains, see Fig. 8. This is done by replacing the reconstructed surface of a given grain with an appropriate FE meshed surface and saving the grains' new geometry 57 Grain-Scale Modeling Approaches for Polycrystalline Aggregates www.intechopen.com with all the imprinted meshes into new ACIS SAT files. The process is repeated for all the grains. Updated ACIS SAT geometry files of individual grains with imprinted surface meshes are imported into a ABAQUS/CAE, assigned appropriate surface definitions, material properties, loads and boundary conditions. Exactly one FE seed per each edge is assigned to preserve the FE meshed surfaces, obtained in the previous step. The number of FE per edge is not allowed to increase/decrease thus automatically creating conformal meshes between adjoining constituents. FE volume-meshing is performed next using ABAQUS/CAE built-in mesher. All the information that has been generated in the previous steps (topology, common surfaces between the constituents, material properties,...) has been saved using the Python pickle module and is now used to hierarchically define all the properties, including node, element and surface sets. Generating FE models of individual grains is independent upon each other which again provides for an efficient parallelization.
In the last step, zero thickness layers of cohesive elements are inserted between the adjacent grains. Layers of zero thickness triangular cohesive elements are inserted between the nodes occupying the same position on the adjacent surfaces. The triangular cohesive elements are oriented to conform with the tetrahedral elements on both surfaces. The nodes, elements, set, surfaces,... and all other definitions are also updated to reflect the new configuration. Fig. 9 illustrates mesh of adjacent grains with inserted zero thickness cohesive layers. Further details on the procedures employed in automatic and parallel generation of the finite element meshes are available in (Simonovski & Cizelj, 2011a). Fig. 10 shows the obtained FE models for the 3D Voronoi tessellations given in Fig. 3. The mesh quality factors are given in Table 1.

Bulk grains
Isotropic elasticity, anisotropic elasticity and anisotropic elasticity with crystal plasticity constitutive laws are commonly used for bulk grains. Since isotropic elasticity can not account for the effects due to different crystallographic orientations of the grains this is not covered here. Overview of the other two constitutive laws is given below.

Anisotropic elasticity
In general, the relation between the elastic stress tensor, σ El. ij , and the elastic strain tensor, ǫ El. kl for anisotropic elasticity is given by Eq. (1).
Here, the C ijkl stands for the stiffness tensor. For materials with a cubic lattice there are only three independent values in the stiffness tensor: C iiii = C 11 , C iijj = C 12 and C ijij = C 44 . The relation between the strains and stresses is then given by Eq. (2).

Crystal plasticity
Crystal plasticity theory (Hill & Rice, 1972;Rice, 1970) assumes that the plastic deformation in monocrystals takes place via a simple shear on a specific set of slip planes. Deformation by other mechanisms such as for example diffusion, twinning and grain boundary sliding is here not taken into account. The combination of a slip plane, denoted by its normal m α i , and a slip direction, s α i , is called a slip system, (α). The plastic velocity gradient,u p i,j , due to a crystallographic slip can be written using Eq. (3) (Needleman, 2000). The summation is performed over all active slip systems,(α), withγ (α) representing the shear rate.
Grain-Scale Modeling Approaches for Polycrystalline Aggregates www.intechopen.com

Grain boundaries with cohesive zone approach
Based on the experimental observations, e.g. (Coffman & Sethna, 2008), grain boundaries are modeled with a cohesive-zone approach using cohesive elements. The traction-separation constitutive behaviour with the damage initiation and evolution as implemented in ABAQUS are used in this work. The cohesive elements are inter-surface elements of often negligible thickness, which essentially measure the relative displacements of the surfaces of adjoining continuum elements. The strains ǫ in the cohesive elements are defined using the constitutive thickness of the element T 0 (mostly different from the geometric thickness which is typically close or equal to zero) and the separations of the element nodes δ as compared to their initial unloaded positions, Eq. (11).
The indices n, s and t denote the normal and two orthogonal shear directions of the cohesive element. The normal direction always points out of the plane of the cohesive element. The tractions on the cohesive elements are then given by Eq. (12).       Typical traction-separation response is given by Fig. 12. Damage evolution D(δ) is defined by Eq. (13) for the normal direction (and both shear directions).
The actual load-carrying capability of the cohesive element in the normal direction would then be [1 − D(δ)] K nn and correspondingly for the two shear directions.

Cohesive elements issues
In this section we explore the response of the cohesive elements in their normal direction. Let us have a cuboid, divided into three grains as depicted in the Fig. 13 Fig. 13. A simple Y model.

64
Polycrystalline Materials -Theoretical and Practical Aspects www.intechopen.com The resulting stress tensor is given by Eq. (19).
The boundaries between the grain are defined with the vectors, normal to the planes of the grain boundaries.
Since we know the stress tensor and the normals for these three planes, we can compute the stresses in the normal direction for each of them. For the plane between the Grain1 and Grain2 Similarly, we obtain the following normal stress values for the other two planes: Grain2Grain3: σ n = 150 MPa (25) Fig. 14 displays the normal stresses in the cohesive elements as calculated from the ABAQUS finite element models. On the left hand side triangular prism cohesive elements are used, whereas on the right hand side rectangular prism cohesive elements are used. One can see that in the case of triangular prisms there is a variation in the normal stress for a given plane, in particular at the Y triple points. For the rectangular prism cohesive elements no such variations are observed and the values from the FE model match exactly with the theoretically computed values. Rectangular prism cohesive elements should therefore be preferentially used.
Unfortunately, meshing complicated shapes with rectangular prism cohesive elements also requires that the grains need to be meshed with rectangular prisms-hexahedral elements.
With shapes as complicated as seen in Fig. 5 this is extremely difficult and one therefore uses triangular prisms-tetrahedral elements. This then requires the use of triangular prism cohesive elements if one is to obtain a conformal mesh between the structural and cohesive elements. One therefore automatically introduces a degree of discrepancy. Fig. 15 displays the computed versus the theoretical normal stresses for the cohesive elements for the 3D Voronoi tessellation with 100 grains with external load of 70 MPa. Isotropic elastic constitutive law is used for the grains with the wire loaded in tension. The solid line represents the ideal response. One can observe a significant scatter from the ideal response. Fig. 16 displays the cohesive elements with red color indicating elements with more than 50 % difference between the theoretical and computed normal stresses. One can see that a significant scatter exists. Similar scatter has been reported on the grain structure simulated using 3D Voronoi tessellations (Kamaya & Itakura, 2009). Vast majority of the problematic cohesive elements are located on the triple lines between the grains. The scatter of the normal stresses of the cohesive elements not lying on the triple lines is significantly smaller, with most values within the ±20 % deviation. Increasing mesh density helps to alleviate the issue to some degree by reducing the area of the problematic elements. Other factors such as the stiffness of the cohesive element and its thickness have negligible effect on the scatter (Simonovski & Cizelj, 2011c).

Examples
In this section we look at the some results for the wire model with initialization and propagation of intergranular stress corrosion cracks. The grain boundaries are modeled using the above described cohesive zone approach and classified into resistant and susceptible grain boundaries, depending upon the crystallographic orientation of the neighboring grains. In this work a simplification is used where resistant grain boundaries are defined as coincidence site lattice (Grimmer et al., 1974) (Σ3 through Σ29) grain boundaries and low angle grain boundaries with misorientation angle between the neighboring grains below 15 o , following (Marrow et al., 2006). Σ value is computed as the ratio enclosed by a unit cell of the coincidence sites and the standard unit cell (Bollman, 1982) from the Rodrigues vectors of the two neighboring grains. Brandon (Brandon, 1996) criterion for a proximity to a coincidence site lattice structure with proportionality constant of 10 o is used. All other grain boundaries are defined as susceptible grain boundaries.
First, we constrain the nodes on the back surface in all three directions. Next, we apply tensile stress of 60 MPa to the wire's front surface. This load is equal to the one in the experiment (King et al., 2008a;. After the application of the load we assume that the wire is exposed to an acidified solution of potassium tetrathionate (K 2 S 4 O 6 ) that penetrates into the wire and degrades the susceptible grain boundaries. This is done to mimic the experiment and is accomplished by employing a user-defined field variable of a cylindrical shape with it's axis aligned with the wire's axis.
The radius of the user-defined field variable cylinder is progressively decreased. Once a cohesive element of a susceptible grain boundary lies outside the user-defined field variable cylinder's radius, its δ 0 n and δδ f n values are decreased, resulting in practically instantaneous full degradation.
Damage of a cohesive element of a susceptible grain boundary that is inside the user-defined field variable cylinder's radius is caused only by mechanical overload, as depicted in Fig. 12. The described approach is purely mechanical. To improve convergence, the decrease of δ 0 n and δδ f n values is not done abruptly. Additionally, a viscous damping of 0.01 is applied to the damage function to improve the convergence during the cohesive elements damage evolution. Susceptible grain boundaries at initial and end 10 % of the wire's length are not allowed to be affected by the user-variable to reduce the edge effect and improve the numerical stability.
The rate at which the radius of the user-defined field variable decreases (i.e. degradation rate) is linked to the stability of the computation. If a large value is selected, then within one computational time increment the separation of the opposite faces in a cohesive element can reach the critical value of δ f n at which the element completely degrades. If this degradation process occurs within one computational time increment, the convergence is degraded, even more so when this occurs simultaneously in several cohesive elements. Small enough degradation rates therefore need to be used, so that the process of degradation of cohesive elements is captured within the computational time increments. Alternatively, one can select very small computational time increments or increase the step time. Similarly, the δ 0 n and δδ f n values should not be very small since at already small load increment the resulting separation of the cohesive element faces could be high enough to instantaneously completely damage the element, again causing convergence issues. Fig. 17 shows the Mises stresses in bulk grains (left part) and damage evolution on the grain boundaries (right part) for the three constitutive models for grains. 30K case with 10.0 µm element size is presented. The first two rows (with the exception of the legend) display the results at an early time increment where the damage due to the corrosion is very limited. For the Mises stress significant differences can be observed between the IE (isotropic elasticity) and AE (anisotropic elasticity) constitutive laws. This is to be expected since the crystallographic orientations are disregarded in IE approach. Due to the low applied load (less than 1/3rd of the yield stress) almost no difference can be observed between the Mises stresses for the AE and AE+CP (anisotropic elasticity+crystal plasticity) constitutive models. This was true even at the end of the simulation where the damage of the grain   boundaries was at its highest, resulting in stress redistribution from the failed areas of grain boundaries to the neighboring areas of grains. The last row therefore displays only the Mises stress for the AE+CP model. Redestribution of the Mises stress can be observed due to the degradation of the susceptible, tensile-loaded grain boundaries which decreases the amount of stress transferred from one side of the boundary to the other. Stress is redistributed in the neighboring areas of grain boundaries and grains, increasing the compressive loading. This can be seen as dark patches in the last row.
The presented approach is, however, not without its deficiencies. First, all tensile-loaded grain boundaries degrade at the same rate. In reality, higher degradation rates might be expected for the susceptible grain boundaries with higher tensile load. Also, the initial grain boundary stiffness is taken to be uniform whereas stiffness distribution based upon the properties of adjacent grains is expected (Coffman & Sethna, 2008). Lastly, the selected corrosive environment penetration approach does not account for the topology of the grain boundary network. These issues will therefore need to be addressed in future work. Slightly better convergence of the AE+CP model was observed. However, the computational times for the AE+CP model were more than twice those for the AE model, see Table 3

Conclusion
With the rapid development of computational capabilities and new experimental techniques we are moving closer to understanding the full role of microstructure on the materials performance. Not only are advanced approaches for simulating microstructures being used but also models of as-measured structures are actively being developed. Among the former, tools such as the presented 2D and 3D Voronoi tessellations can be used whereas for the latter, experimental techniques such as the X-ray diffraction contrast tomography which enable 3D characterization of grains are indispensable. Basics of these approaches are covered here. Both approaches share a common difficult task of creating a finite element model in terms of both appropriate meshes and model sizes. Surface reconstruction issues and complex geometry make the process more difficult. The presented approach effectively deals with some of these issues. The others, however, remain and are subject to further work and research.
The demonstration of the approach is presented on several cases of 3D Voronoi tessellations and two cases of a 400 µm diameter stainless steel wire. In all cases the grain boundaries are explicitly modeled using the cohesive zone approach with zero physical thickness finite elements. Grain boundaries are classified into resistant and susceptible grain boundaries, depending upon the crystallographic orientation of the neighboring grains. Grain boundary damage initialization and early development is then computed for a stainless steel case for several constitutive laws, ranging from isotropic elasticity up to crystal plasticity for the bulk grain material. Since isotropic elasticity approach disregards the crystallographic orientations it should not be used in these cases. Little differences were observed between the anisotropic elasticity and anisotropic elasticity+crystal plasticity approaches, with the latter resulting in more than twice as long computation times.
In all cases almost uniform degradation of the grain boundaries is observed. This is attributed to a) a missing link between the grain boundary load and rate at which the corrosion penetrates the grain boundaries, b) uniform grain boundary stiffnesses whereas stiffness distribution based upon the properties of adjacent grains is expected (Coffman & Sethna, 2008) and c) the selected corrosive environment penetration approach does not account for the topology of the grain boundary network. These issues are the subject of further work.
The numerical stability of the simulation including damage is reasonable, with slightly better convergence for the anisotropic elasticity+crystal plasticity approach. The degradation of a cohesive element is linked to the stability of a simulation. If this degradation process occurs within one computational time increment, the convergence is degraded, even more so when this occurs simultaneously in several cohesive elements. Small computational time increments should therefore be used. Similarly, the δ 0 n and δδ f n values should not be very small since at already small load increment the resulting separation of the cohesive element faces could be high enough to instantaneously completely damage the element, again causing convergence issues.