Automated Reconstruction and Conforming Mesh Generation for Polycrystalline Microstructures from Imaging Data

Abstract

A new algorithm named PolyCISAMR is introduced to automatically generate high-fidelity conforming finite element (FE) meshes for two-dimensional polycrystalline microstructures. PolyCISAMR extends the capabilities of the Conforming to Interface Structured Adaptive Mesh Refinement (CISAMR) algorithm, which transforms a structured grid overlaid on the domain geometry into a high-quality conforming mesh. The PolyCISAMR approach uses a segregated meshing strategy, where CISAMR is used to discretize each grain independently and the resulting matching meshes are merged to form the final FE model. In addition, this article presents a set of integrated algorithms for processing low-resolution images of a polycrystal, reconstructed using DREAM.3D software (Version 6.5.121), to generate NURBS characterizations for each grain prior to mesh generation. Example problems demonstrate the effectiveness of PolyCISAMR in creating high-quality meshes for various polycrystalline metallic microstructures along with corresponding crystal plasticity finite element (CPFE) simulations.

1. Introduction

Polycrystalline materials are composed of several grains with sizes that range from several nanometers to millimeters and typically have random crystallographic orientations [1]. Most inorganic materials, including all common metals, most ceramics, rocks, and ice, are polycrystalline in nature and are composed of a large number of crystallites joined together by grain boundaries [2]. The mechanical behavior of polycrystalline materials is highly dependent on their microstructure, specifically the size and orientation of grains, as well as the area of grain boundaries. For instance, the Hall–Petch relationship states that the strength of polycrystals increases as the average grain size decreases [3]. Additionally, the manufacturing process plays a crucial role in determining the microstructural features of polycrystalline materials. This includes factors such as the cooling rate and the extent of deformation, which can greatly impact macroscopic properties of the material such as strength and toughness. Various numerical techniques have been developed to study the micromechanical response of polycrystals, including mean-field approaches (e.g., Taylor and self-consistent schemes) [4,5,6], Crystal Plasticity Finite Element (CPFE) [7], and Green’s function Fast Fourier Transform (FFT) models [8,9]. While CPFE and FFT models are more computationally demanding than mean-field approaches, they can accurately capture the evolution of individual grain shapes, the interaction between neighboring grains, intra-granular strain evolution, and the heterogeneous onset of damage [10].

CPFE models provide the advantage of incorporating complex microstructural features and heterogeneities, such as micro-cracks, inclusions, and precipitates, as well as enabling simulations with complex boundary conditions resulting from inter/intra-granular interactions [11,12]. However, creating image-based finite element (FE) models can be challenging. It is crucial to accurately represent grain morphologies, including grain boundaries and corners, when creating a high-fidelity FE model for CPFE simulations. Low resolution or the presence of noise in pixelated/voxelated imaging data obtained from sources such as electron backscatter diffraction (EBSD) can greatly hinder the achievement of this goal. Virtually reconstructed microstructures using methods such as Voronoi tessellation [13,14] or DREAM.3D [15] can lack accuracy or lead to similar challenges when generating high-fidelity FE models. Reconstruction algorithms relying on the Voronoi tessellation method approximate the grain boundaries using straight planar surfaces, which is not realistic for many complex polycrystalline microstructures. Additionally, DREAM.3D synthesizes grain models with jagged boundaries (similar to EBSD data), and does not capture their edge/corner geometry. To bypass these challenges, CPFE simulations typically use pixelated/voxelated meshes [13,14]. However, these meshes cannot accurately represent stress concentrations at the grain boundaries, which can negatively impact the accuracy of simulations aimed at modeling damage/debonding along grain boundaries.

Several meshing algorithms are implemented to generate high-quality conforming FE meshes, such as Delaunay triangulation [16,17], Quadtree/Octree-based approaches [18,19], and advancing front methods [18,20,21]. The Delaunay triangulation method solves a nonlinear system of equations to optimize the mesh node locations, with the objective of ensuring that no mesh node is located within the circumcircle of a triangular element in 2D or the circumference of a tetrahedral element in 3D. In the advancing front method, new triangular/tetrahedral elements are created on an active front that emanates from specific regions in the domain, usually the boundaries, and gradually advance towards the interior to discretize the whole domain. Quadtree/Octree-based methods recursively refine quadrilateral/hexahedral elements to reach a desired level of refinement, resulting in conforming elements with arbitrary aspect ratios along material interfaces and domain boundaries. Mesh smoothing or optimization-based approaches are then employed to improve the element aspect ratios by iteratively relocating mesh nodes. While these methods can generate high-quality meshes for complex geometries/microstructures, they encounter two main challenges: (i) the high computational cost of iterative techniques such as mesh smoothing that are used to improve the element aspect ratios [22]; and (ii) occasional difficulty in convergence for intricate geometries, which might necessitate user intervention during the meshing process to change meshing parameters or manipulate the domain geometry.

There have been several attempts at adapting or customizing meshing algorithms to generate conforming meshes for polycrystals. Neper [23,24] is an open-source polycrystal reconstruction and meshing package that uses Voronoi/Laguerre tessellation techniques to mesh synthesized microstructures. Other researchers have also adopted Voronoi/Laguerre tessellation-based approaches for meshing polycrystals, as seen in [25,26]. Although these iterative techniques can generate conforming meshes for simplified polycrystalline microstructures (e.g., with planar grain boundaries), their implementation for modeling more complex and realistic polycrystals can be challenging. DREAM.3D [15] can create more sophisticated synthetic polycrystalline microstructures with realistic grain shapes, size distribution, and crystallographic orientations compared to Neper. However, due to the challenges associated with meshing such complex polycrystals, DREAM.3D is only equipped with the ability to generate voxelated meshes. Although suitable for CPFE-based analyses, the use of such meshes can result in unrealistic stress concentrations and lower accuracy due to jagged representation of the grain boundaries. Pyle et al. [27] compared CPFE simulations on voxelated and conforming meshes and found that while macroscopic behavior, texture evolution, and stress/slip statistics may be similar, the conforming grain boundary representation provides a more accurate misorientation distribution and smoother stress field. For modeling intergranular damage, such as fatigue crack growth, generating a conforming mesh is crucial for simulating the problem. An open-source package worth mentioning for creating grain-conforming meshes is the XtalMesh Toolkit [28], which is capable of generating high-fidelity FE models for realistic 3D polycrystalline microstructures reconstructed using DREAM.3D.

This study presents a comprehensive computational framework for transforming pixelated polycrystalline microstructure images into high-fidelity 2D conforming meshes. Before generating the mesh, a robust preprocessing engine is required to transform the raw pixelated image (here, synthesized using DREAM.3D) into an appropriate geometrical model of the polycrystal with smooth representation of the grain boundaries. This approach utilizes advanced image processing and computational geometry algorithms to identify individual grains, boundaries, and vertices, and characterize grain morphologies using Non-Uniform Rational B-Splines (NURBS). A modified version of the Conforming to Interface Structured Adaptive Refinement (CISAMR) algorithm [29,30,31], coined PolyCISAMR, is developed to generate conforming meshes for the resulting NURBS-based model. The CISAMR is a non-iterative mesh generation algorithm that transforms a structured mesh overlaid with the domain geometry into a conforming mesh. The proposed PolyCISAMR algorithm combines the non-iterative meshing capability of CISAMR with the added advantage of independent meshing of each grain, making it a segregated meshing algorithm highly suitable for parallel implementation. The novelty of the proposed research is two-fold: (1) we introduce an image processing algorithm to extract the NURBS representation of grains from pixelated imaging data, providing smooth representation of grain boundaries and characterizing their sharp corners; and (2) a segregated version of the CISAMR algorithm is introduced for generating conforming meshes for the resulting NURBS-based model of a polycrystal, in which each grain is discretized independently.

The remainder of this manuscript is structured as follows. In Section 2, we present the preprocessing algorithm used to convert the raw image of a polycrystal to a NURBS-based geometrical model. The PolyCISAMR algorithm and its implementation details are introduced in Section 3. The governing equations for crystal plasticity are presented in Section 4, where we describe the CPFE model as well. In Section 5, several example problems are presented to demonstrate the ability of PolyCISAMR to generate high-quality conforming meshes and corresponding CPFE simulations, including a mesh convergence study and a comparison with pixelated meshes. Finally, the conclusions of this work are summarized in Section 6.

2. Reconstruction of Polycrystalline Model

2.1. Synthesizing Polycrystyals Using DREAM.3D

DREAM.3D (Version 6.5.121) is an open-source software environment that is used for microstructure quantification and representation across various material classes and length scales [15]. It is an effective tool for the virtual reconstruction of polycrystalline microstructures with realistic grain shapes, size distributions, and crystallographic orientations. The software allows users to export synthesized polycrystals in a variety of formats, including two-dimensional Portable Network Graphics (PNG) image slices, Stereolithography (STL) files containing triangulated surface meshes for individual grains, and eXtensible Data Model and Format (XDMF) files. However, the PNG images generated by DREAM.3D do not have sufficient resolution for creating realistic grain morphologies and identifying sharp edges/corners. Similarly, even after using the Laplacian smoothing filter, the STL files of individual grains generated by the software have a jagged representation of grain boundaries and cannot be used to create accurate FE models or approximate stress concentrations during CPFE simulations. To overcome these limitations, we use the XDMF format to export microstructures from DREAM.3D into the software package ParaView [32] for visualization. As shown in Figure 1, 2D image slices can then be exported from ParaView in the form of high-resolution PNG images (in this case, with a resolution of $1100\mathrm{pixels}\times 1100\mathrm{pixels}$). Although this approach still produces a jagged representation of grain boundaries, it has sufficient resolution for further image processing and creation of a smooth boundary representation using NURBS patches.

2.2. Identifying Individual Grains

An algorithm utilizing color-based boundary detection was developed to identify individual grains in 2D image slices of synthesized polycrystalline microstructures, as shown in Figure 2. The algorithm starts by applying a color-based mask to the image in order to separate grains of similar color. A pixel erosion algorithm is then implemented to identify grain boundaries. This process involves recursively visiting each grain, eroding it by a few pixels in all directions, and then taking its dot product with the white background after binarization. This process is illustrated in Figure 2 for one of the grains belonging to the polycrystalline microstructure depicted in Figure 1. Next, a skeletonization algorithm [33] is employed to extract a pixel-wide skeleton representation of the grain boundaries, as shown in Figure 3. After successfully identifying each grain, the algorithm proceeds to determine its neighboring grains, which is crucial in order for later stages of the preprocessing phase to identify shared grain boundaries. This is achieved by creating an enlarged bounding box (BBox) around each grain and using it to quickly identify neighboring grains that overlap with it, meaning that at least one of their pixels falls inside the BBox.

The presence of extremely small features (grains) composed of a limited number of pixels in 2D image slices can lead to difficulties during mesh generation, including the need for excessive mesh refinement and the formation of degenerate elements. Therefore, prior to identifying grain boundaries, it is crucial to identify and eliminate any small, unrealistic features composed of only a few pixels in order to avoid such challenges. To address this issue, the first step is to use an area threshold criterion to identify small grains. Principal Component Analysis (PCA) [34] is then applied to detect any small grains that exhibit a skewed orientation. As seen in Figure 4, such elongated small grains can negatively impact the meshing process by causing unacceptably high element aspect ratios if not treated properly beforehand. Unlike non-elongated small grains, a higher threshold must be used for the area calculation to eliminate elongated small grains. It is important to note that in PCA, the point cloud of a small grain serves as input to determine its principal components. A large difference between principal components indicates asymmetry or possible elongation of a grain. Excessively small and small, elongated grains are eliminated by merging them with one of their neighboring grains, chosen in a way that results in higher convexity after the collapse, as shown in Figure 4. It is important to note that the convexity is evaluated as the ratio of the region’s area to its convex hull [35]. This measure ranges from 0 to 1, where 1 refers to a perfectly convex polygon.

2.3. NURBS Characterization of Grains

To improve the rough appearance of grain boundaries in the polycrystal image and create a suitable input file for the PolyCISAMR mesh generation algorithm, NURBS functions are used to describe the geometry of each grain. A NURBS curve $\mathbf{C}\left(u\right)$ is a rational generalization of B-splines $M_i^p\left(u\right)$ which are functions of a parametric coordinate u. A knot vector $\mathbf{U}=u_1,u_2,u_3,\dots ,u_{n+p+1}$ represents the parametric coordinates of the initial and end points of a NURBS patch in a $C^0$-continuous NURBS curve. In this knot vector, n and p refer to the number of control points and the order of B-Splines, respectively. A NURBS curve can be written as [36]

$$\mathbf{C}\left(u\right)=\sum _{i=1}^n\frac{{\mathbf{x}}_i{w}_i{M}_i^p\left(u\right)}{{\sum }_{j=1}^n{w}_j{M}_j^p\left(u\right)},$$

(1)

where ${\mathbf{x}}_i$ is a set of control points (described in the Cartesian coordinate system) and $w_i$ are weights of each control point.

To evaluate the NURBS characterization of each grain, we must first identify its sharp corners. A junction detection algorithm [37] is implemented that receives the skeleton image $I_{\mathrm{skel}}$ of the polycrystal (cf. Figure 3) as the input and detects the sharp corners, as shown in Figure 5. These sharp corners are in fact intersection points of two or more grain boundaries in a given 2D image slice. The junction detection process begins by applying a 2D convolution operation [38] on the binary skeleton image $I_{\mathrm{skel}}$ (cf. Figure 5a) using a $3\times 3$ box blur kernel of all ones [39]. The resulting image $I_{\mathrm{blur}}$, which is depicted in Figure 5b, has pixel values equal to the sum of pixel values in the Moore neighborhood of each pixel of $I_{\mathrm{skel}}$ [40]. A pixel-wise dot-product is then taken between $I_{\mathrm{skel}}$ and $I_{\mathrm{blur}}$, which generates an image with possible pixel values of $\{0,3,4,5,6,7,8,9\}$ depending on the number of grain boundaries that intersect at a given pixel, as shown in Figure 5c. Here, a value of 0 denotes the interior of a grain, 3 represents pixels on the grain boundaries, and any value greater than 3 corresponds to pixels near the grain vertices (sharp corners). Next, pixel values greater than 3 are filtered out from $I_{\mathrm{blur}}$ to obtain multiple junctions containing sharp corners. These junctions can either have a pixel value greater than 3 (junction B in Figure 5c) or multiple pixels that meet this threshold (junction A). A sharp corner in the latter scenario is identified by calculating the mean of all pixel locations and then selecting the closest pixel to this mean value (Figure 5d).

Next, vertices that are too close to domain boundaries or extremely close to each other are detected; these are respectively snapped to the boundary or collapsed to their midpoint (see Figure 6). While identifying such vertices based on their distance to other sharp corners is straightforward, there are cases where multiple sharp corners can be in close proximity, as shown in Figure 7. Such scenarios require a special approach to ensure that small grains are preserved in the final polycrystal model after merging the sharp corners. Note that after generating the skeleton image of the initial polycrystal we no longer have access to individual grains at this stage of processing, as only grain boundary information is available and it is unfeasible to identify small grains that must be preserved. Therefore, as depicted in Figure 7, for each sharp corner we first determine the number of neighboring sharp corners $N_{\mathrm{s}}$ within a distance $l_{\mathrm{gs}}$ given by

$$l_{\mathrm{gs}}=\frac{l+w}{2N},$$

(2)

where l and w are the length and width of the 2D polycrystal image and N is the total number of grains. If $N_{\mathrm{s}}>3$, it implies that this sharp corner belongs to a small grain relative to the average size of grains in the polycrystalline microstructure and should not be relocated/merged.

To obtain the NURBS characterization of grain boundaries, we first implement the Shi–Tomasi corner detection algorithm [41] to identify a set of corner points that represent grain boundaries. Note that the detection of corner points here refers to points that result from the jagged representation of the grain boundaries (cf. Figure 7), and should not be confused with the sharp corners detected via the junction detection algorithm discussed previously. To avoid any ambiguity, we refer to such corner points (detected on a jagged boundary) as feature points in the remainder of this section. The Shi–Tomasi corner detection algorithm [41] starts by identifying small image windows/patches that have large variations in intensity (large gradients). The algorithm computes the change in intensity $E(u,v)$ of each image patch $I(x,y)$ when moved slightly in the x and y directions by u and v, respectively, i.e.,

$$E(u,v)=\sum _{x,y}^{}{[I(x+u,y+v)-I(x,y)]}^{2}w(x,y),$$

(3)

where $w(x,y)$ is a weighting function. In the simplest case, $w(x,y)=1$; however, it could be a Gaussian-like function instead to emphasize the central area of the image patch. Applying the Taylor series expansion to (3) yields

$$E(u,v)\approx [u,v]\mathbf{M}\left[\begin{array}{c}u\\ v\end{array}\right],$$

(4)

where $\mathbf{M}$ is a symmetric $2\times 2$ matrix of the system provided by

$$\mathbf{M}=\sum _{x,y}^{}w(x,y)\left[\begin{array}{cc}{I}_x{I}_x& I_x{I}_y\\ I_x{I}_y& I_y{I}_y\end{array}\right].$$

(5)

In Equation (5), $I_x$ and $I_y$ represent the image gradients in the x and y directions, respectively. A feature point on the jagged boundary is identified by a large variation of $E(u,v)$ in all directions of $[u,v]$. Because the eigenvalues of M represent the variance in the corresponding principal directions, this matrix should have two large eigenvalues, ${\lambda }_1$ and ${\lambda }_2$, at a feature point. However, because the minimum of these two eigenvalues represents the smallest amount of variance in the image patch, it is less affected by noise and is a more reliable criterion for feature point detection [41]. Therefore, a scoring function R is defined as $R=min({\lambda }_1,{\lambda }_2)$. If this scoring function is greater than a threshold ${\lambda }_{\mathrm{th}}$, it is determined that the image window contains a feature point. Note that this threshold ${\lambda }_{\mathrm{th}}$ is chosen as a fraction of the best score $R_{\max }$ from a given image, i.e., ${\lambda }_{\mathrm{th}}=q{R}_{\max }$, where we have chosen $q=0.01$.

After identifying the feature points along grain boundaries (smaller green points in Figure 8), a subset of these points is selected as the NURBS control points, which is shown as larger blue points in Figure 8. Note that using all feature points as control points leads to spurious oscillations in the resulting NURBS representation of grain boundaries. To select appropriate control points, we loop over the global set of feature points, then for each point being visited we find the neighboring points within the distance $l_{gs}$, as defined in (2). All identified points are then removed, and only the original point being visited is maintained as a control point. Figure 9 illustrates the control points and corresponding NURBS curves resulting from this algorithm, representing the pixelated polycrystalline microstructure depicted in Figure 1. However, as described next, the construction of these NURBS curves is not simple and requires further processing, as individual grains and their boundaries are unidentified at this stage and we only have two sets of points representing the entire polycrystalline microstructure (sharp corners of grains and control points on their boundaries).

To identify individual grains, we first use the Moore-neighbor contour tracing algorithm with Jacob’s stopping criterion [42] to obtain a sorted set of points along grain boundaries (cyan points in Figure 10). This algorithm begins at one of the black boundary pixels of a grain boundary (the start pixel) in the skeleton image and explores its Moore neighborhood (eight adjacent pixels) in the counterclockwise direction [40] to find the next boundary pixel. This algorithm sorts a set of pixels that represent an individual grain counterclockwise, with Jacob’s stopping criterion ensuring that the algorithm terminates when the tracer re-enters the start pixel. After identifying/sorting pixels of all grains, we must identify the NURBS control points belonging to each grain. First, the control points located inside the bounding box (BBox) of each grain, enlarged by a factor of $1.5l_{\mathrm{gs}}$, are identified (dark blue and yellow points in Figure 10). Next, we identify control points within the distance $l_{\mathrm{gs}}$ from the sorted grain boundary points previously obtained using the Moore-neighbor contour tracing algorithm. This step is essential to filter out control points belonging to the current (visited) grain (dark blue points in Figure 10) and exclude those of adjacent grains (yellow points). Sharp corners of each grain are identified similarly.

The next step is to identify highly convex grains using a convexity measure $k_{\mathrm{conv}}$ [43], which is the ratio of the grain’s area to the area of its convex hull (the smallest convex polygon containing all of its sharp corners and control points [35]). For grains that do not deviate very much from their convex hull, identified with $k_{\mathrm{conv}}>0.85$, the sharp corners and control points are sorted in a counterclockwise order based on the angle they make with the grain’s centroid. Sorting these points allows the grain geometry to be characterized as a set of NURBS curves between its consecutive sharp corners, as shown in Figure 9. However, when $k_{\mathrm{conv}}<0.85$ (cf. Figure 11), using a centroidal angle-based algorithm to sort the sharp corners and control points may yield erroneous results due to the concave shape of the grain. In such scenarios, we again employ the boundary points obtained using the Moore-neighbor contour tracing algorithm, which are already sorted in counterclockwise order (the cyan points in Figure 10), to sort the control points. Figure 12 illustrates the final NURBS representation of the polycrystalline microstructure image depicted in Figure 3. The pseudocode of the algorithm for reconstructing this NURBS model described above is summarized in Algorithm 1.

Algorithm 1 (Preprocessing phase for image-to-NURBS conversion in PolyCISAMR)
function image_to_nurbs_conversion($X_{\mathsf{3}\mathsf{D}}$)
I ← get_2D_image ($X_{\mathsf{3}\mathsf{D}}$)	▹ get 2D image slices (I) from 3D polycrystalline microstructure ($X_{3\mathrm{D}}$)
grain_set, $I_{\mathsf{GB}}$ ← identify_individual_grains (I)	▹ identify individual grains and grain boundaries ($I_{\mathrm{GB}}$) (cf. Figure 2)
$I_{\mathsf{skel}}$ ← skeletonize ($I_{\mathsf{GB}}$)	▹ get a skeleton ($I_{\mathrm{skel}}$) representation of grain boundaries (cf. Figure 3)
sharpcorners ← detect_junctions ($I_{\mathsf{skel}}$)	▹ get sharpcorners/ grain vertices of all grains (cf. Figure 5)
sharpcorners ← preserve_small_grains (sharpcorners, $l_{\mathsf{gs}}$)	▹ identify small grains and preserve them (cf. Figure 7)
feature_points ← get_feature_pts ($I_{\mathsf{skel}}$)	▹ get feature points for all grain boundaries (cf. green points in Figure 8)
control_points ← get_control_points (feature_points, $l_{\mathsf{gs}}$)  ▹ get control points by removing redundant feature points within a radius of $l_{\mathrm{gs}}$
for ($i=1$ to $\mathsf{size}\mathsf{(}$grain_set)) do	▹ loop over all grains
$B_i$ ← trace_boundary (grain${}_i$)	▹ get a sorted set of points ($B_i$) representing the grain boundary of current grain
BBox${}_i$ ← get_enlarged_BBox (grain${}_i$, $1.5l_{\mathsf{gs}}$)	▹ get an enlarged BBox around the current grain
neighbor_grains${}_i$ ← get_neighbor_grains (grain_set, BBox${}_i$)	▹ get neighboring grains of the current grain
$P_i$ ← filter_control_points (BBox${}_i$, control_points, $l_{\mathsf{gs}}$) ▹ Filter control points ($P_i$) belonging to current grain (cf. Figure 10)
$S_i$ ← filter_sharp_corners (BBox${}_i$, sharpcorners, $l_{\mathsf{gs}}$) ▹ Filter sharp corners/vertices ($S_i$) belonging to current grain (cf. Figure 10)
${\chi }_i$ ← calculate_convexity (grain${}_i$)	▹ Calculate convexity (${\chi }_i$) of the current grain
if (${\chi }_i>0.85$) then
$C_i$ ← calculate_grain_centroid (grain${}_i$)	▹ find centroid ($C_i$) of the current grain
$S_i$ ← sort_sharp_corners ($S_i$, $C_i$)	▹ sort sharp corners counter-clockwise based on centroidal angle
$P_i$ ← sort_control_points ($P_i$, $C_i$)	▹ sort control points counter-clockwise based on centroidal angle
else
$S_i$ ← sort_sharp_corners ($S_i$, $B_i$)	▹ sort sharp corners based on boundary points
$P_i$ ← sort_control_points ($P_i$, $B_i$)	▹ sort control points based on boundary points
end if
$f_i$ ← write_NURBS_file (neighbor_grains${}_i$, $S_i$, $P_i$)	▹ write .nurbs file for the current grain
end for
end function

3. PolyCISAMR Algorithm

3.1. Overview of CISAMR Algorithm

The PolyCISAMR algorithm introduced here for meshing polycrystalline microstructure is an extension of the 2D CISAMR algorithm originally introduced in [44]. CISAMR is a non-iterative algorithm that transforms an initial structured mesh composed of rectangular elements into a high-quality conforming mesh consisting of rectangular and triangular elements in three stages: h-adaptive refinement of elements near material interfaces, r-adaptivity of nodes of elements that intersect the interface, and sub-triangulation of elements with hanging nodes and elements deformed during the r-adaptivity phase. The process is illustrated in Figure 13, and each step of the algorithm is described in more detail in the remainder of this subsection.

∘

h-adaptive refinement: Using NURBS curves as an explicit representation of material interfaces, we first identify four-node quadrilateral (Q4) elements of the background mesh intersecting each interface. As shown in Figure 13a, a Structured Adaptive Mesh Refinement (SAMR) algorithm is then implemented to refine these elements as well as their selected neighboring elements. The objectives of this step are twofold: (i) to reduce the geometric discretization error; and (ii) to reduce the error associated with the gradient recovery in the vicinity of the interface in the corresponding FE simulation. At each refinement level, all nonconforming Q4 elements that are cut by the interface are subdivided into four Q4 sub-elements. The neighboring elements that are subjected to refinement are selected based on the relative distances of their nodes to the interface. In this approach, if the interface intersects an element edge with length h at a distance $d<0.5h$ from one of the nodes of that edge, all neighboring elements that share this node are refined. As discussed in [44], a number of constraints must be satisfied during the SAMR phase to ensure that the final conforming elements have proper aspect ratios. For example, as depicted in Figure 13a, each edge of a background element is allowed to hold only one hanging node. If more hanging nodes appear on an edge, the elements sharing that edge are subjected to additional refinement until this constraint is satisfied. Similarly, adaptive refinement is activated when the edges of an element intersect either an interface multiple times or with multiple interfaces.

∘

r-adaptivity: The next step is to visit the nodes of adaptively refined elements cut by an interface to determine whether they must be relocated to that interface. Figure 13b shows the deformed mesh created after applying this r-adaptivity phase, during which ≈50% of the originally nonconforming elements of the background mesh are transformed to conforming elements. We assume that the lengths of edges connected to mesh node $\mathcal{N}$ are $h_i$ and the intersection point(s) of these edges with the interface have distance(s) $d_i$ from $\mathcal{N}$. As shown in Figure 14, the new location of this node after performing r-adaptivity is determined as follows:

• Case 1: If only one of the edges connected to $\mathcal{N}$ intersects an interface:

  (a)

  If $d\ge 0.5h$: no need to relocate the node.

  (b)

  If $d<0.5h$: move the node to the edge-interface intersection point.

• Case 2: If two of the edges connected to $\mathcal{N}$ intersect an interface:

  (a)

  If $d_1\ge 0.5h_1$ and $d_2\ge 0.5h_2$: do not relocate the node.

  (b)

  If $d_1<0.5h_1$ and $d_2\ge 0.5h_2$: move the node to the closer intersection point at distance $d_1$ and ignore the presence of the intersection point at distance $d_2$.

  (c)

  If $d_1<0.5h_1$ and $d_2<0.5h_2$ and $d_1<d_2$: move the node to the closer intersection point at distance $d_1$ and discard the farther intersection point at distance $d_2$.

∘

Sub-triangulation: Finally, all background elements that were deformed during the r-adaptivity process, along with elements with hanging nodes created due to h-adaptive refinement of their neighboring elements, are subdivided into conforming three-node triangular (T3) elements. The final conforming mesh created after applying this step is depicted in Figure 13c. Note that all remaining nonconforming elements after the completion of the r-adaptivity phase are diagonally cut by the interface. To minimize the aspect ratio of the resulting T3 sub-elements, the following rules are applied while subdividing Q4 elements (cf. Figure 15):

• Case 1: If a Q4 element does not intersect the interface along the diagonal emanating from its smallest angle ${\theta }_{min}$, T3 sub-triangles are created by cutting that element along the diagonal corresponding to its largest angle.
• Case 2: Otherwise, there are two possible scenarios:
  • If ${\theta }_{min}>{60}^{\circ }$: subdivide the element by cutting along ${\theta }_{min}$ to ensure that the aspect ratios of the resulting T3 elements are acceptable.
  • If ${\theta }_{min}<{60}^{\circ }$: subdivide the element by cutting along both diagonals, resulting in the creation of four conforming sub-triangles.

It is worth mentioning that the CISAMR algorithm guarantees an upper bound of three on-element aspect ratios for problems with smooth interfaces, which might increase to four for domain boundaries and interfaces with sharp corners (e.g., the polycrystalline microstructures discussed next). More algorithmic details are provided in [44], along with the pseudocode of the original CISAMR algorithm for creating hybrid Q4-T3 meshes.

3.2. Meshing an Individual Grain Using CISAMR

In the PolyCISAMR algorithm, the CISAMR technique is first employed to discretize individual grains as independent meshing problems, followed by merging the resulting meshes to build the final FE model. However, due to the presence of grains’ sharp corners, creating a high-quality conforming mesh for each grain requires adding new algorithmic features to CISAMR. Because the sharp corners of each grain ${\mathcal{C}}_i$ are already identified during the preprocessing phase (see Section 2.3), we first implement a BBox-based quad-tree search to quickly locate the background elements that contain each sharp corner. In most cases where the grain does not have an acute angle at the sharp corner, a small modification to the CISAMR algorithm is sufficient to generate an appropriate conforming mesh. As illustrated in Figure 16, after performing h-adaptivity (SAMR), the closest mesh node to the sharp corner is identified and snapped to that. The subsequent phases of the CISAMR algorithm (r-adaptivity and sub-triangulation) are carried out in the same way as described in Section 3.1 to construct the conforming mesh.

For sharp corners corresponding to an acute angle, which can be categorized into three different scenarios, the following additional treatments are required (cf. Figure 17):

• Case 1: As shown for the highlighted background element in Figure 17a, although this element intersects the grain, all of its nodes are located outside the grain boundaries. This would be a special case, as the relative location of a background element with respect to the grain boundaries (inside, outside, on the boundary) is determined based on the relative location of its nodes. Therefore, it is important to ensure that such cases are not misclassified during CISAMR implementation, as any element for which all four nodes are outside the grain is typically categorized as an outside element. Fortunately, because such special cases always occur in the vicinity of the background element containing the sharp corner, it is easy to identify such misclassified elements. In this approach, if all four nodes of the background element containing a sharp corner fall outside the grain, we perform an additional check on the relative location of its neighboring elements to ensure that at least one of them intersects the grain boundaries. If this is not the case (similar to the special case discussed here), the location of midpoints, quarter points, etc., of all edges of these elements relative to the grain boundaries is recursively checked until an intersecting element is identified.
• Case 2: This scenario, which is also depicted in Figure 17a and often accompanies Case 1, occurs during the r-adaptivity phase. Here, the background element edge identified with nodes ${\mathcal{N}}_1$–${\mathcal{N}}_2$ intersects with the grain boundary twice and both intersection points are in the same half of the edge, i.e., $d<0.5h$ for both intersections points relative to ${\mathcal{N}}_1$. Therefore, because node ${\mathcal{N}}_2$ does not relocate and node ${\mathcal{N}}_1$ can only relocate to one of these intersection points (the closest intersection point), the resulting grain geometry would be unrealistic without any treatment (see Figure 17a). The treatment is rather simple: introducing a new midpoint node (hanging node) on edge ${\mathcal{N}}_1$–${\mathcal{N}}_2$, which can then be relocated to the other intersection point during r-adaptivity. As shown in the “after treatment” configuration in Figure 17a, this simple modification yields an appropriate conforming mesh while accurately capturing the grain geometry.
• Case 3: Similar to Case 2, in this scenario a background mesh node must be relocated to two intersection points, this time on two different edges connected to that and not on the same edge, as shown for ${\mathcal{N}}_1$ in Figure 17b). In such cases, picking either of the intersection points to relocate the mesh nodes leads to misrepresentation of the grain geometry, as shown in the without treatment scenario in Figure 17b. To handle such cases, a new hanging node is added on the edge with the intersection point that has a larger distance to ${\mathcal{N}}_1$ (here, edge ${\mathcal{N}}_1$–${\mathcal{N}}_2$) before performing r-adaptivity. Relocating this hanging node to the intersection point followed by performing the sub-triangulation phase of CISAMR allows the grain geometry to be accurately captured in the final mesh.

Another necessary modification in the CISAMR algorithm for meshing individual grains is the treatment of background elements with a very small edge in the vicinity of sharp corners after performing r-adaptivity, as shown for edge ${\mathcal{N}}_1$–${\mathcal{N}}_2$ in Figure 18. Without proper treatment, the sub-triangulation of such excessively deformed background elements leads to the creation of sub-triangles with an unacceptably high aspect ratio. Such cases can be identified by checking the ratio of the minimum element edge length $h_{min}$ connected to a sharp corner (after r-adaptivity) to the original background mesh size $h_0$. If $h_{min}/h_0<\xi$, where $\xi =0.25/2^k$ and k is the number of SAMR levels, we apply an edge collapse algorithm to that edge by snapping the node that is not located on the sharp corner to the corner node. In Figure 18, this treatment strategy leads to collapsing ${\mathcal{N}}_2$ to ${\mathcal{N}}_1$, which eliminates the small background element edge and results in elements with proper aspect ratios after applying the sub-triangulation phase.

3.3. PolyCISAMR Segregated Meshing Algorithm

The PolyCISAMR algorithm relies on a segregated meshing strategy, meaning that it uses CISAMR to independently mesh individual grains without any information exchange (communication between processors in the parallel implementation) during the meshing process. The resulting meshes, which are automatically assured to have conforming nodes along their edges, are then stitched together by merging these nodes to build a conforming mesh for the polycrystalline microstructure. The main steps of this algorithm, which are schematically illustrated in Figure 19, are as follows:

• Step 1: After superimposing the polycrystalline microstructure onto a structured background mesh, a subset of this mesh that overlaps with the slightly enlarged BBox of each grain is selected. The BBox is enlarged by $0.1h$ in each direction, where h is the size of the background elements. It is worth noting that in the parallel implementation (and even for the sequential version), generating a large background mesh and extracting subsets of it for each grain would not be not necessary; instead, a small structured mesh that overlaps with the enlarged BBox of each grain can directly be created, provided that the nodal coordinates of this mesh snap to nodes of an imaginary structured grid covering the entire domain.
• Step 2: As described in Section 3.2, the CISAMR is then implemented to create a conforming mesh for each grain as an independent meshing problem using the background mesh identified in the previous step. The only restriction at this step is that the same level of SAMR must be used for meshing grains in order to achieve conforming meshing along the edges of neighboring grains. In the parallel implementation, each grain is assigned to a different partition and meshed independently, with no communication between processors during this process. This distinct feature is the reason for referring to PolyCISAMR as a “segregated meshing algorithm”. Note that because all operations that lead to creating new nodes (SAMR), node relocation (r-adaptivity), and sub-triangulation in these segregated meshes are non-iterative and identical during the CISAMR implementation, it is automatically guaranteed that the resulting meshes will have matching nodes along adjacent grain boundaries.
• Step 3: The segregated meshes generated in Step 2 for different grains are then stitched together (merged) to construct the final FE model. In this process, we iterate over all nodes of the individual grain meshes and store them in a red–black tree-based associative container [45], with the coordinates as key and the node pointer as value. In this approach, we do not need to identify neighboring grains and sort the nodes located on their boundaries in a specific order during the merging process. Instead, we can easily search for and identify matching nodes along the grain boundaries through a simple coordinate check and then merge them by assigning the same global node ID in the final mesh. Note that the tree-based associative container has a time complexity of $O\left(log\right(n\left)\right)$ for insertion/search operations, and does not allow duplicate keys (node coordinates), allowing this task to be performed efficiently.

The PolyCISAMR algorithm described here has several advantages, both in its sequential version, as presented in this work, and when implemented in parallel. The ability to mesh individual grains one by one as an independent problem allows for the handling of massive microstructures composed of hundreds or even thousands of grains without facing memory restrictions. Figure 20 illustrates the conforming mesh generated using this algorithm for the NURBS representation of the polycrystalline microstructure depicted in Figure 12. The background meshes generated for each grain correspond to a $120\times 120$ background grid for the entire domain with a grid size $h=1.67\mathsf{\mu }$m. Using two levels of SAMR along each grain boundary, the final mesh consists of 63,432 elements and 43,113 nodes, with a maximum element aspect ratio of 3.8. Note that only two elements in this mesh have an aspect ratio larger than three due to the sharp angle at one of the grain corners. The pseudocode for the PolyCISAMR algorithm is presented in Algorithm 2.

Algorithm 2 (PolyCISAMR algorithm)
function Generate_Single_Grain_Meshes($f_{\mathsf{list}}$, SAMR_levels)
for ($i=1$ to $\mathsf{size}\mathsf{(}{f}_{\mathsf{list}}$)) do	▹ loop over all the nurbs files ($f_{\mathrm{list}}$)
Bounds${}_i$, $S_i$, $P_i$ ← read_nurbs_file ($f_i$)	▹ get bounds, sharp corners ($S_i$), control points ($P_i$) for the current grain
BBox${}_i$ ← get_enlarged_BBox (Bounds${}_i$, $0.1h$)	▹ get an enlarged Bounding Box around current grain
back_mesh${}_i$ ← get_background_mesh (BBox${}_i$)	▹ get background mesh for the current grain
$M_i$ ← SAMR (back_mesh${}_i$, $P_i$, SAMR_levels)	▹ perform h-adaptive refinement
$M_i$ ← r_adaptivity_sharp ($S_i$, $M_i$)	▹ locate nearest nodes to the sharp corners ($S_i$) and snap them (cf. Figure 16)
$M_i$ ← capture_sharp_interfaces ($M_i$) ▹ Additional treatment for extremely sharp interfaces during r-adaptivity (cf. Figure 17)
$M_i$ ← regular_r_adaptivity ($M_i$)	▹ Perform regular r-adaptivity
$M_i$ ← collapse_small_edges ($M_i$, $\xi$)	▹ Collapse small edges with edge length ratio less than $\xi$ (cf. Figure 18)
$M_i$ ← sub_triangulation ($M_i$)	▹ Perform sub-triangulation (cf. Figure 15)
end for
end function
function Merge_Single_Grain_Meshes($M_{\mathsf{list}}$)
red_black_tree<coords,node_pointer> node_tree	▹ Container storing unique point and node pointer
for ($i=1$ to $\mathsf{size}($$M_{\mathsf{list}}$)) do	▹ Loop over grain sub-meshes ($M_{\mathrm{list}}$)
for ($j=1$ to $\mathsf{size}\left(\mathsf{el}\_\mathsf{set}\right)$) do	▹ Loop over each element ($e_j$) in a sub-mesh
if is_outside$\left(e_j\right)$ then	▹ Filter elements not lying inside a grain
continue
end if
for ($k=1$ to $\mathsf{size}\left(\mathsf{nodes}\right)$) do	▹ Loop over each node ($n_k$) in an element $e_j$
if is_in_node_tree$(n_k$.coords())==false then	▹ If mesh point is not in $\mathrm{node}\_\mathrm{tree}$
node_tree←add_to_node_tree$(n_k.\mathsf{coords}\left(\right),n_k)$	▹ Add to $\mathrm{node}\_\mathrm{tree}$
$N_{\mathsf{glob}}\leftarrow$ assign_global_nodeID$\left(n_k\right)$	▹ Assign node ID for final merged mesh
end if
$N_{\mathsf{glob}}\leftarrow$ get_global_nodeID$(n_k.\mathsf{coords}\left(\right))$	▹ Get node ID for a mesh point already in the $\mathrm{node}\_\mathrm{tree}$
end for
$E_{\mathsf{con}}\leftarrow$ store_element_connectivity($N_{\mathsf{glob}}$) ▹ Store element connectivity for the current element using global node IDs
end for
end for
write_nodes$\left(N_{\mathsf{glob}}\right)$	▹ Write nodal data in mesh input file
write_elem_connectivity$\left(E_{\mathsf{con}}\right)$	▹ Write element connectivity in mesh input file
end function

4. Crystal Plasticity Model

In this section, we review the general framework of the kinematics and constitutive equations of crystal plasticity, as presented in [46,47]. Additionally, we specialize the multi-mechanism model for Grade 91 steel and the two-phase material system of Zircaloy-2 with hydrides, which were originally presented in [48,49], respectively.

Kinematic Relations

The mesoscale behavior of polycrystalline metallic alloys can be well-described by the CPFE method [50,51,52], which postulates that plastic deformation occurs through the sliding of dislocations on preferential slip systems with respect to the crystal lattice orientation in each material grain. To represent the anisotropic plastic strain and the potentially large lattice rotation, various kinematic treatments of the constitutive equations have been proposed using hyperelastic (strain energy-based) [51,53] or hypoelastic (objective stress rate-based) [50,54] approaches. In this case, the Green–Naghdi hypoelastic approach is adopted from [47] due to its advantage of formulating the stress update formulas in the unrotated frame [55], which separates the large strain kinematics from the shear stress to slip rate relations, allowing for the incorporation of multiple models in a unified implementation. The Green–Naghdi objective rate of the Cauchy stress tensor $\mathit{\sigma }$ is expressed through the material time derivative $\left(\dot{•}\right)={\displaystyle \frac{\partial (•)}{\partial t}}$ as

$$\check{\mathit{\sigma }}=\dot{\mathit{\sigma }}+\mathit{\sigma }\Omega -\Omega \mathit{\sigma }.$$

(6)

The right polar decomposition of the deformation gradient $\mathbf{F}=\mathbf{RU}$ into rotation tensor $\mathbf{R}$ and stretch tensor $\mathbf{U}$ is used to define the skew-symmetric angular velocity tensor $■=\dot{\mathbf{R}}{\mathbf{R}}^T$. For crystal plasticity of metals and minerals with stiff elastic modulus, the elastic strains $\mathit{\epsilon }\ll \mathbf{I}$ are quite small compared to the identity tensor $\mathbf{I}$. Therefore, the traditional multiplicative split of $\mathbf{F}$ into elastic and plastic parts is assumed to take a simplified form through polar decomposition of each component and retention of the small elastic strains, i.e.,

$$\mathbf{F}={\mathbf{F}}^{\mathrm{e}}{\mathbf{F}}^{\mathrm{p}}\approx \left(\mathbf{I}+\mathit{\epsilon }\right){\mathbf{R}}^{\mathrm{e}}{\mathbf{R}}^{\mathrm{p}}{\mathbf{U}}^{\mathrm{p}}.$$

(7)

This treatment accommodates features of plastic anisotropy on slip systems and finite elastic rotation ${\mathbf{R}}^e$. The derivations in [47] provide the substitution of (7) into (6) and mapping to the unrotated intermediate configuration to yield the rate equation for the unrotated Cauchy stress $\mathbf{t}={\mathbf{R}}^T\mathit{\sigma }\mathbf{R}$ provided by

$$\dot{\mathbf{t}}={\mathbf{C}}_0:\left(\mathbf{d}-{\mathbf{d}}^{\left(\mathrm{p}\right)}\right)+{\mathbf{Rw}}^{\mathrm{p}}{\mathbf{R}}^T\mathbf{t}-{\mathbf{tRw}}^{\mathrm{p}}{\mathbf{R}}^T.$$

(8)

The first term in this equation involves the constant anisotropic elastic fourth-order elasticity tensor ${\mathbf{C}}_0$ in the lattice frame multiplying the elastic strain rate $\dot{\mathit{ϵ}}$. This is additively decomposed into the total unrotated rate of deformation $\mathbf{d}=\frac{1}{2}\left(\dot{\mathbf{U}}{\mathbf{U}}^{-\mathbf{1}}+{\mathbf{U}}^{-\mathbf{1}}\dot{\mathbf{U}}\right)$ and the plastic deformation rate ${\overline{\mathbf{d}}}^{\mathrm{p}}$. The other two terms involving the plastic vorticity ${\overline{\mathbf{w}}}^{\mathrm{p}}$ are correcting terms due to plastic anisotropy. The major postulate of crystal plasticity theory is that plastic straining occurs due to the accumulation of slip of dislocations on preferred crystallographic planes. Hence, the plasticity tensors ${\overline{\mathbf{d}}}^{\mathrm{p}}$ and ${\overline{\mathbf{w}}}^{\mathrm{p}}$ are expressed via plastic slip rates ${\dot{\gamma }}^{\left(s\right)}$ resolved onto primary slip planes $s=1,\dots ,n_{\mathrm{slip}}$

$$\begin{array}{cc}\hfill {\overline{\mathbf{d}}}^{\mathrm{p}}& =\sum _{s=1}^{n_{\mathrm{slip}}}{\dot{\gamma }}^{\left(s\right)}{\mathbf{R}}^{\mathrm{p}T}{\tilde{\mathbf{m}}}^{\left(s\right)}{\mathbf{R}}^{\mathrm{p}},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\overline{\mathbf{w}}}^{\mathrm{p}}& =\sum _{s=1}^{n_{\mathrm{slip}}}{\dot{\gamma }}^{\left(s\right)}{\mathbf{R}}^{\mathrm{p}T}{\tilde{\mathbf{q}}}^{\left(s\right)}{\mathbf{R}}^{\mathrm{p}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{m}}}^{\left(s\right)}& =\mathrm{sym}\left({\tilde{\mathbf{b}}}^{\left(s\right)}\otimes {\tilde{\mathbf{n}}}^{\left(s\right)}\right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{q}}}^{\left(s\right)}& =\mathrm{skew}\left({\tilde{\mathbf{b}}}^{\left(s\right)}\otimes {\tilde{\mathbf{n}}}^{\left(s\right)}\right),\hfill \end{array}$$

(12)

where ${\tilde{\mathbf{b}}}^{\left(s\right)}$ denotes the slip direction unit vector and ${\tilde{\mathbf{n}}}^{\left(s\right)}$ is the slip plane unit normal. The curved “tilde” overbar indicates that they are in the lattice (reference) frame. The slip rates ${\dot{\gamma }}^{\left(s\right)}$ depend upon the (unrotated) Cauchy stress and a set of $n_{\mathrm{hard}}$ internal hardening variables $\mathit{\xi }={\xi }^{\left(1\right)},{\xi }^{\left(2\right)},\dots ,{\xi }^{\left(n\right)}$, written generically as

$${\dot{\gamma }}^{\left(s\right)}={\dot{\gamma }}^{\left(s\right)}(\mathbf{t},\mathit{\xi }),\dot{\mathit{\xi }}=\dot{\mathit{\xi }}(\mathbf{t},\mathit{\xi }).$$

(13)

Meanwhile, the plastic rotation evolves via the plastic vorticity tensor ${\overline{\mathbf{w}}}^{\mathrm{p}}$ according to

$$\dot{{\mathbf{R}}^{\mathrm{p}}}={\overline{\mathbf{w}}}^{\mathrm{p}}{\mathbf{R}}^{\mathrm{p}}.$$

(14)

At each material point, the coupled set of time-continuous Equations (8) and (13) govern the evolution of the primary unknowns $\mathbf{t}$ and $\mathit{\xi }$. To model them numerically, we adopt the implementation from [46] that employs backwards Euler time discretization for all quantities except the plastic rotation (14) and solves the coupled system with Newton–Raphson supplemented by line search techniques. A salient feature of this approach is that it is modular with respect to the slip rate and hardening models through the linkages (13), for which two specific models are provided in the following sections. The reader is encouraged to consult [46,48] for linearization and time discretization details.

5. Numerical Examples

In this section, the quality of the mesh generation algorithm is evaluated through a series of CPFEM simulations performed on a Grade 91 ferritic–martensitic steel and a Zircaloy-2 alloy (constitutive models are summarized in Appendix A and Appendix B). The polycrystalline microstructure was reconstructed virtually using DREAM.3D, assuming an average grain size of 33 $\mathsf{\mu }$m [48], a random crystal texture, $\mu =3.5$, $\sigma =0.1$, and a primary equiaxed shape. Multiple two-dimensional image slices were extracted from the resulting microstructure and transformed into NURBS models using the image processing algorithms presented in Section 2. The PolyCISAMR algorithm introduced in Section 3 was then employed to generate the conforming meshes in the form of input files for the Warp3D CPFE code [56]. Because this code only performs 3D analyses, the meshes were first extruded into a single layer of quadratic wedge elements and subdivided into three quadratic tetrahedral elements for the CPFE analysis.

In the following simulations, the boundary conditions on the CPFEM models are consistent: symmetry is applied to the left, lower, and back faces of the unit cell, traction-free conditions are applied to the top y and front z faces, and a constant axial true strain rate of 0.01 1/hr is applied to the right x face. This strain rate is larger than the slower creep strain rates of 0.0001 1/hr for which the model was calibrated; however, the macroscale stress–strain curves produced from the model at this rate are on a similar order of magnitude as the yield stress of 300 to 400 MPa for Gr91 at this temperature and strain rate [48]. The time step size varies from 0.02 to 0.05 hr during the analysis.

5.1. Mesh Refinement Study

First, we investigate the accuracy of stress calculations on three meshes generated using PolyCISAMR with varying levels of refinement on background meshes with $h=7.5$ $\mathsf{\mu }$m (coarse mesh), $5.0$ $\mathsf{\mu }$m (medium mesh), and 2.5 $\mathsf{\mu }$m (fine mesh) on one of the image slices (Slice 1). The medium mesh is illustrated in Figure 21. All meshes are generated using two levels of SAMR, resulting in a total number of 46,500, 75,721, and 189,448 elements, respectively. In none of these meshes does the maximum aspect ratio of the elements exceed 3. More details regarding the characteristics and computational cost associated with the construction of these meshes are provided in

Table 1. Total mesh generation time ($t_{\mathrm{tot}}$), individual grain meshing time ($t_{\mathrm{g}}$), and merging time ($t_{\mathrm{m}}$) using PolyCISAMR for meshing Slice 1 model with different refinement levels, where $N_{\mathrm{n}}$ and $N_{\mathrm{e}}$ denote the number of nodes and the number of elements, respectively.

Mesh Resolution	${\mathit{N}}_{\mathbf{n}}$	${\mathit{N}}_{\mathbf{e}}$	${\mathit{t}}_{\mathbf{tot}}$	${\mathit{t}}_{\mathbf{g}}$	${\mathit{t}}_{\mathbf{m}}$	${\mathit{r}}_1={\mathit{t}}_{\mathbf{tot}}/{\mathit{N}}_{\mathbf{e}}$	${\mathit{r}}_2={\mathit{t}}_{\mathbf{s}}/{\mathit{N}}_{\mathbf{e}}$
Coarse	34,044	51,534	22.55 s	21.02 s	1.53 s	$4.38\times {10}^{-4}$	$4.08\times {10}^{-4}$
Medium	56,910	84,382	32.53 s	29.72 s	2.81 s	$3.86\times {10}^{-4}$	$3.52\times {10}^{-4}$
Fine	148,413	206,024	74.86 s	65.31 s	9.55 s	$3.63\times {10}^{-4}$	$3.17\times {10}^{-4}$

. The time required for transforming the initial pixelated image to a NURBS model (input file for mesh generation) was 15 s. The total PolyCISAMR mesh generation time ($t_{\mathrm{tot}}$) scales in a superlinear fashion with the number of elements ($N_{\mathrm{e}}$) for meshes with different refinement levels (see $r_1=t_{\mathrm{tot}}/N_{\mathrm{e}}$ values for different meshes). Similarly, according to Table 1, the time required for constructing individual grain meshes ($t_{\mathrm{g}}$) scales superlinearly with $N_{\mathrm{e}}$. This is attributed to the relatively small time required for merging segregated grain meshes ($t_{\mathrm{m}}$) compared to the overall computational time $t_{\mathrm{tot}}$, as reported in Table 1. Additionally, this table shows that $t_{\mathrm{m}}$ scales well with the number of nodes ($N_{\mathrm{n}}$) in the mesh. These performance metrics demonstrate the efficiency of the PolyCISAMR algorithm for meshing large polycrystalline domains composed of many grains.

Plane stress CPFE simulations under uniaxial tension were conducted on these meshes. The resulting stress, total strain, plastic strain, and hardening fields for the fine mesh are illustrated in Figure 22 and Figure 23. Two applied strain levels were considered: a lower value of 0.0026, at which some grains begin to develop plastic flow, and a larger value of 0.02, at which most grains have hardened and reached saturation stress. At the first load level, there is significant variation in the axial stress ${\sigma }_{xx}$ contour shown in Figure 22a, as certain grains are favorably oriented for plastic slip while others are not. On the other hand, at the higher applied strain level in Figure 23, the axial stress is more uniform across the domain and mostly equal to the saturation stress of 310 MPa. It is important to note that these contour plots are not projected nodal stress averages; they are the elemental stress averaged over the quadrature points. The lack of color discontinuities in each figure indicates the high quality of the mesh.

Analyzing the results for all three meshes showed that the homogenized (volume-averaged) stress–strain curves were nearly identical, indicating that even the coarsest mesh is sufficiently refined to accurately capture the macroscale response. However, local quantities are expected to converge more slowly, which was indeed the observation made in this study. For example, Figure 24 shows the effect of mesh refinement levels on the axial stress in the subset of Slice 1.

5.2. Comparison with Pixelated Meshes

In this section, we compare the two CPFEM simulations performed on PolyCISAMR conforming meshes with their pixelated counterparts. As noted previously, pixelated meshes are widely used in the literature for CPFEM studies due to the easy process of creating such meshes [13,14]. However, the jagged or stair-stepped representation of grain boundaries in pixelated meshes leads to the introduction of artificial stress concentrations in these regions, which can undermine the accuracy of CPFEM simulations [27]. In the comparison study presented in this section, CPFEM results on two PolyCISAMR meshes generated on a background mesh with $h=5$ $\mathsf{\mu }$m, one with no SAMR (18,705 elements) and another with 1 level of SAMR (36,680 elements), are compared to simulations conducted on pixelated meshes with $h=5$ $\mathsf{\mu }$m (15,625 elements) and $h=2.5$ $\mathsf{\mu }$m (62,500 elements), respectively. Small portions of these meshes are illustrated in Figure 25, which clearly shows the difference in the representation of grain boundaries in conforming versus pixelated meshes. Note that, unlike in the conforming meshes, the geometrical representation of grains in the pixelated meshes depends on the mesh size; although the mesh refinement reduces the stair-stepped appearance of grain boundaries, it does not completely eliminate this artifact, as shown in Figure 26.

The CPFEM approximations of the variation of normal stress ${\sigma }_{xx}$ in the PolyCISAMR conforming and pixelated meshes at an applied strain of 0.02 are illustrated in Figure 25. As shown in the inset of each figure, these plots show element-wise values of ${\sigma }_{xx}$ in each mesh, which converge more slowly compared to a smooth stress field obtained using nodal averaging or more sophisticated techniques. Therefore, rather than comparing the stress values, studying the pattern of stress variations would be more useful in understanding the effect of mesh type (conforming vs. pixelated) on the accuracy. The drastic difference between simulated normal stresses in the PolyCISAMR and pixelated meshes can clearly be seen in the highlighted regions marked A, B, and C in the inset of the meshes shown in Figure 25. Note the unrealistic stress concentrations in the pixelated meshes compared to the conforming meshes in these regions, which emanate from the jagged representation of grain boundaries in the former. Although a different window of the microstructure is inspected in Figure 24, even on the coarse conforming PolyCISAMR mesh in Figure 24a there are no oscillations in color such as are prevalent in Figure 25. The comparison between the highlighted region C in the conforming mesh with one level of SAMR (cf. Figure 25b) and the fine pixelated mesh (cf. Figure 25d) is especially interesting, as it shows a single element in the latter with erroneously high stress due to the inaccurate representation of the grain boundary geometry in this region. It is worth noting that such artifacts of a pixelated mesh only affect the stress/strain approximation in a small portion of elements of the whole mesh along grain boundaries; thus, their impact on averaged values such as homogenized stress-strain curves is muted. However, such localized errors could be devastating when approximating phenomena such as intergranular damage using cohesive elements, which requires accurate recovery of stresses along grain boundaries.

5.3. Effect of Random 2D Slices

Next, we examine the impact of the polycrystalline microstructure on the computed mechanical fields when different 2D slices are selected from the same 3D microstructure. First, five 2D slices were randomly selected from a 3D microstructure reconstructed using DREAM.3D. Note that while the statistics of different 2D slices may not be identical to the 3D model, they are relatively similar. The grain orientations in each slice were chosen at random and the number of grains in each case may vary slightly. Figure 27a,c illustrates the microstructure and the conforming meshes generated using PolyCISAMR for three of the five slices analyzed in this section, which consist of 1.8 to 2 million elements. The maximum aspect ratio of elements in all five meshes is lower than 3.9.

First, we present the macroscale stress–strain curve for various slices in Figure 28. All of the slices exhibit similar trends for the initial yield point and the shape of the curve. As shown in the inset of this figure, even at the saturation point, the curves are still very similar. Next, Figure 29 illustrates the axial stress contour plot at a 0.02 strain for three slices. Although the contour values at a specific point $(x,y,z)$ vary greatly between slices, the overall distribution of low and high stresses appears similar, with roughly two to four grains having a low (blue) value and another two to four having a high (red) value. Additionally, stress gradients are mostly found along grain boundaries, with only two to four grains displaying diagonal bands of color gradients on the interior.

5.4. Effect of Including Hydride Precipitates

In the last example, we investigate the impact of including hydrides in the polycrystalline microstructure on the mesh quality and CPFEM results. Zircaloy is selected as the material system due to the presence of needle-like hydrides formed during the manufacturing process as well as in-service irradiation. To accurately represent Zircaloy, a synthetic microstructure was created in DREAM.3D using the parameters $\mu =4.0$ and $\sigma =0.2$ and a Primary Equiaxed shape. A 2D slice with dimensions of $950\mathsf{\mu }\mathrm{m}\times 950\mathsf{\mu }\mathrm{m}$ was then extracted, converted to a NURBS model, and meshed using the PolyCISAMR algorithm. The resulting conforming mesh, which consists of $\approx 0.25$ million elements with a maximum aspect ratio of 3.2, is illustrated in Figure 30. The precipitates introduced along grain boundaries represent 5.4% of the volume fraction. A random texture was sampled for grain orientations, with a hexagonal close-packed (HCP) crystal structure. The hydride phase was modeled using an isotropic plasticity model with higher stiffness and strength compared to the Zircaloy matrix. The microstructural model was subjected to plane stress tensile loading, as in the previous examples, although the strain rate was lower at 0.000025 1/hr to match the calibrated conditions from [49] for the constitutive model (see Appendix B).

The CFPE approximation of the normal stress and strain fields in this slice at two different applied strains are illustrated in Figure 31 and Figure 32. Compared to previous examples, a notable increase in sites of stress concentrations in the alloy microstructure is evident in this problem, especially at an applied strain of 0.008. This behavior is attributed to the presence of the hydride phase, which is stiffer and attracts more stress due to deforming less than the surrounding matrix. Additionally, the HCP Zircaloy phase is more anisotropic than the BCC Grade 91 material from the previous studies in both the elastic and plastic regimes. This greater anisotropy leads to larger stress gradients between grains with different crystal orientations, which is another contributing factor to the localizations. Despite the presence of sharp gradients in the fields, a close examination of the contours indicates that the fields remain free of oscillations and are fairly well-resolved on the high-quality conforming meshes generated using PolyCISAMR.

6. Conclusions

The CISAMR algorithm has been extended to generate conforming meshes for polycrystalline microstructures in 2D. This non-iterative algorithm, named PolyCISAMR, uses a segregated meshing strategy to discretize individual grains independently of other grains and then merge them to build the FE model, which makes it suitable for meshing huge microstructures and allows for efficient parallelization in the future. An automated image processing algorithm has been presented for converting low-resolution pixelated images of a polycrystal to a NURBS representation of grains for use as an input file to the PolyCISAMR meshing algorithm. This manuscript further addresses new algorithmic aspects of CISAMR such as the treatment of extremely sharp corners of grains by using an element collapse approach to preserve the mesh quality. A detailed mesh refinement study and several example problems have been presented to demonstrate the application of the proposed framework in modeling polycrystalline microstructures. Two alloy systems were adopted from the previous modeling work of the co-authors, and the computed macroscale trends were in agreement with these prior results. Additionally, in stark contrast to the pixelated meshes that we investigated, the elemental contour fields computed from the high-quality meshes generated using PolyCISAMR were free of oscillations and captured stress localizations near grain boundaries. Key insights from this work include:

• Our numerical examples show that the presented framework is suitable for modeling polycrystalline microstructures.
• The CPFE simulations performed on conforming meshes generated using PolyCISAMR are superior to those obtained using pixelated meshes, as the artificial stress concentrations along the jagged boundaries in the latter are eliminated when using PolyCISAMR conforming meshes.

Future work involves the extension of the proposed image processing and segregated meshing algorithms to 3D, which requires further modifications to the latest version of the CISAMR algorithm [57] to ensure conformity of the nodes along neighboring grain boundaries, especially along their sharp edges/corners. When this is accomplished, it will be possible to merge individual grain meshes to generate conforming meshes for the 3D polycrystals. Moreover, as mentioned earlier, the PolyCISAMR algorithm in its current form is segregated in nature, enabling efficient parallelization without any communication overhead. Future efforts will be devoted to developing a parallel version of PolyCISAMR, which can then be used to mesh massive polycrystalline microstructures.