A novel 3D anisotropic Voronoi microstructure generator with an advanced spatial discretization scheme

The input of the algorithm consists of (i) an initial hexahedral mesh, covering the domain of interest Ω, (ii) a set of N_b domain boundary functions ${f\hspace{-2pt}}_{1}(\overrightarrow{y}\enspace ),\dots ,{f}_{{N}_{\mathrm{b}}}(\overrightarrow{y}\enspace )$ and their gradients $\overrightarrow{\nabla }{f\hspace{-2pt}}_{1}(\overrightarrow{y}\enspace ),\dots ,\overrightarrow{\nabla }{f}_{{N}_{\mathrm{b}}}(\overrightarrow{y}\enspace )$ whose definition will be provided later, (iii) the seed point coordinates ${\overrightarrow{x}}_{1},\dots ,{\overrightarrow{x}}_{N}$, (iv) the corresponding preferred growth direction vectors ${\overrightarrow{p}}_{1},\dots ,{\overrightarrow{p}}_{N}$ (including optional second preferred growth vectors ${\overrightarrow{q}}_{1},\dots ,{\overrightarrow{q}}_{N}$) and (v) the corresponding volumes of the ellipsoids V₁, ..., V_N . The initial hexahedral mesh can be a regular voxel mesh, but this is not required; it can also be an irregular hexahedral mesh which has been generated with standard hexahedral finite element discretization algorithms in case the domain of interest is more complex. Together with the initial hexahedral mesh, domain boundary functions ${f\hspace{-2pt}}_{k}(\overrightarrow{y}\enspace )$ should be provided. Each of these functions ${f\hspace{-2pt}}_{k}(\overrightarrow{y}\enspace )$ are defined to be zero at the kth face of the domain, Γ_k , reading

$$\begin{equation}\begin{cases}{f\hspace{-2pt}}_{k}(\overrightarrow{y}\enspace )=0\quad & \quad \forall \enspace \overrightarrow{y}\in {{\Gamma}}_{k},\\ {f\hspace{-2pt}}_{k}(\overrightarrow{y}\enspace )\ne 0\quad & \quad \forall \enspace \overrightarrow{y}\notin {{\Gamma}}_{k}.\end{cases}\end{equation} \tag{ 8 }$$

For example, signed, unsigned or squared distance functions fulfill these conditions. The use of these boundary functions is twofold. First, for each node p in the mesh with position ${\overrightarrow{y}}_{p}$, let K_p contain the set of domain boundary face numbers k on which this node is located, i.e.

$$\begin{equation}{K}_{p}=\left\{k,\vert f{\hspace{-2pt}}_{k}({\stackrel{\to }{y}}_{p})\vert < {\epsilon}\right\},\end{equation} \tag{ 9 }$$

with a given tolerance, and define B_p = #K_p as the number of boundaries in the set K_p . The boundary functions can then be used to classify nodes in the mesh based on their location, i.e. either in the interior (B_p = 0), on a domain face (B_p = 1), on a domain edge (B_p = 2) or at a domain vertex (B_p = 3). Second, the boundary functions are used to displace nodes on particular faces while ensuring that these nodes remain on the corresponding faces. The remaining input consists of data (${\overrightarrow{x}}_{i}$, ${\overrightarrow{p}}_{i}$, ${\overrightarrow{q}}_{i}$ and V_i ) representing the characteristics of each grain. This can be based on experimental observations or statistics of the (morphology of the) microstructure. In order to represent a specific part of a microstructure, it is also possible to extract the required input from (3D) electron backscatter diffraction (EBSD) data by fitting ellipsoids to grains. One can also obtain the preferred growth direction vector ${\overrightarrow{e}}_{p}$ as the largest temperature gradient from e.g. a thermal process simulation, since crystals tend to grow in that direction.

For the elaboration of the detailed numerical implementation, the example shown in figure 2 is considered. The initial hexahedral mesh consists of a regular voxel mesh of 10 × 10 × 5 elements having a total size of $1\times 1\times \frac{1}{2}$. The corresponding six signed distance boundary functions (one for each face of the domain of interest depicted in figure 2) are given as

$$\begin{equation}\begin{aligned}{f\hspace{-2pt}}_{1}(\overrightarrow{y}\enspace )& =\overrightarrow{y}\cdot {\overrightarrow{e}}_{1},\\ {f\hspace{-2pt}}_{2}(\overrightarrow{y}\enspace )& =\overrightarrow{y}\cdot {\overrightarrow{e}}_{2},\\ {f\hspace{-2pt}}_{3}(\overrightarrow{y}\enspace )& =\overrightarrow{y}\cdot {\overrightarrow{e}}_{3},\\ {f\hspace{-2pt}}_{4}(\overrightarrow{y}\enspace )& =1-\overrightarrow{y}\cdot {\overrightarrow{e}}_{1},\\ {f\hspace{-2pt}}_{5}(\overrightarrow{y}\enspace )& =1-\overrightarrow{y}\cdot {\overrightarrow{e}}_{2},\\ {f\hspace{-2pt}}_{6}(\overrightarrow{y}\enspace )& =\frac{1}{2}-\overrightarrow{y}\cdot {\overrightarrow{e}}_{3}.\end{aligned}\end{equation} \tag{ 10 }$$

In section 3.3, this choice is motivated in more detail. Furthermore, the five considered grain seed points along with their ellipsoidal velocity growth fields have been depicted in figure 2 as well.

Zoom In Zoom Out Reset image size

The first part of the algorithm is to assign each element to a specific grain. This is done in two steps. First, for each element e the time required to grow from each grain seed point i to that element is determined using equation (2), in which $\overrightarrow{y}={\overrightarrow{y}}_{e}$ is the center of element e. Based on these growth times ${t}_{1}({\overrightarrow{y}}_{e}),\dots ,{t}_{N}({\overrightarrow{y}}_{e})$, the grain seed points are sorted from minimum to maximum growth time. Preferably, each element e is eventually assigned to the first grain seed point in this list. For the example case, the preferred assignment of elements to grains is shown in figure 3(a), in which the color of the element corresponds to the seed point velocity field in figure 2. Furthermore, the light-shaded color of the elements corresponds to the preferred grain seed point assignment (which is the result of the current, first step) and a dark-shaded element corresponds to the actually assigned grain seed point (which will be detailed in the second step). By considering the red and purple grain in figure 3(a), it is clear how the anisotropic, ellipsoidal growth speed of the grains is taken into account. In this first step, interactions between grains are not taken into account. Consequently, a grain may be split into multiple, unconnected parts by another grain. In the example case (figure 3(a)), the yellow grain consists of two unconnected parts (it is also not connected in the interior). However, the grain could not have grown from the seed point on the right to the left part, since growth would be blocked by the purple grain in between. Because the grain seed point is on the right, the left part will be referred to as an 'island'. The left part is thus physically not representative and unacceptable. Therefore, the preferred assignment of grains cannot be accepted directly.

Zoom In Zoom Out Reset image size

In the second step of the assignment of elements to grains, it is ensured that two important aspects are not present in the final grain connectivity: grain islands and non-manifold meshes. The former has been discussed in the previous paragraph and the importance of the latter will be discussed in section 3.3. An element e is identified as part of an island when there is no path of neighboring elements assigned to the same grain that connects e to the element containing the corresponding seed point. Here, and in the rest of the detailing of the algorithm, elements are only defined as neighbors when they share a common face. This is achieved by iteratively accepting the preferred assignment of grains. During the first iteration, the elements that contain a grain seed point are accepted if their preferred assignment corresponds to the same grain seed point. This means that in case two or more seed points are located within the same element, only the grain seed point that grows to the center of that element first, is accepted. In figure 3(b), the state of the example case is shown after the first iteration. Again, the light-shaded color of the elements corresponds to the color of the preferred grain seed point assignment and a dark-shaded element corresponds to the actually assigned grain seed point. During all subsequent iterations, the preferred assignment of element e to grain i is accepted when (i) at least one neighboring element has already been accepted and assigned to the same grain i as element e and (ii) it does not introduce a non-manifold condition in the connectivity of grain i. A non-manifold condition is obtained when the connectivity of a given mesh contains a feature that cannot exist in the 'real' world, e.g. a position where elements are only connected via a single vertex or edge. In order to identify non-manifold grain connectivity, the Euler number E, defined as

$$\begin{equation}E={N}_{\text{nodes}}-{N}_{\text{edges}}+{N}_{\text{faces}},\end{equation} \tag{ 11 }$$

is utilized. Given a 3D hexahedral mesh with N_nodes boundary nodes, N_edges boundary edges and N_faces boundary faces, the Euler number E must be equal to 2 in order to be free of non-manifold conditions [22]. For each node p of element e, the Euler number ${E}_{p,i}^{e}$ is calculated for the set of elements that share node p and are assigned to grain i, including element e. Furthermore, the Euler number ${E}_{p,{\neg}i}^{e}$ (considering the set of elements that share node p and are not assigned to grain i, thus excluding element e) is also computed. If for any p it holds that ${E}_{p,i}^{e}\ne 2$ or ${E}_{p,{\neg}i}^{e}\ne 2$, a non-manifold condition would be introduced upon accepting the preferred assignment of grain i to element e, which is therefore rejected. In figures 3(c)–(e), the state of the example case is visualized for iterations 2–4, respectively.

After a number of iterations, it may occur that no new elements are accepted in the last iteration. However, it is possible that there are still elements that have not been accepted yet. This is the case in figure 3(f), where the state of the example case is depicted after 18 iterations. In the next iterations, elements are also allowed to be assigned to their second preferred grain seed point, which, for the example case, have been indicated in figure 3(g). A number of iterations later, only a single element, as depicted in figure 3(h), has not yet been accepted by the algorithm. This is because accepting this grain seed point assignment would introduce a non-manifold condition in the connectivity of the blue grain (it is not an island since it can be reached from below the element, see figure 3(i) in which only the elements corresponding to the blue grain are depicted). Therefore, in the final grain connectivity shown in figure 3(j), this element has been assigned to the closest grain seed point that does not introduce an island or non-manifold condition, which is the purple grain. Using this algorithm, the final grain connectivity does not contain of any grain islands or non-manifold conditions.

For specific applications that rely on element regularity, e.g. using an fast Fourier transform (FFT) solver, the resulting voxel-based mesh can be used directly. However, in many other applications, e.g. when using crystal plasticity and a finite element solver, a non-grain-conforming mesh is undesirable as it can introduce spurious stress concentrations [23]. Therefore, the next step is to convert the mesh to a grain-conforming discretization. This is accomplished using four steps, which are detailed in the following.

First, so-called pillow elements are inserted at the grain boundary interface as proposed by Owen et al [22]. Therefore, faces shared by two elements that have been assigned to two different grain seed points are extruded to both sides of the interface. Each grain boundary face thus introduces two new hexahedral elements, effectively surrounding each grain by a single layer of new elements, hence the name pillow elements. In figure 4(a), the resulting mesh is shown for the example case. A new node p^⋆ corresponding to grain i is obtained by (duplicating and) displacing its parent node on the grain boundary p by a distance d in the local grain boundary inward normal direction ${\overrightarrow{n}}_{p}$. This normal vector can be obtained by averaging the inward normals of the grain boundary element faces that share node p. Notice that for this step it is crucial to have a mesh free of non-manifolds. By using the local grain boundary normal direction (and not for example the direction vector pointing to the geometric center of the grain), concave shaped grains can also be handled by the algorithm, as e.g. illustrated by the green or blue grain in figure 4(a). The pillow size d can be chosen arbitrarily (small), however, it must be smaller than half the local mesh size. Otherwise, the nodes of an element with pillowing on opposite faces would overlap, which directly introduces an invalid element. For visualization purposes, in figure 4(a), d = 0.02, which is $\frac{1}{5}$ of the mesh size. The main purpose of this step is to ensure that each element has maximum one face on a grain boundary interface. Consequently, the grain conformity is automatically addressed in the element connectivity.

Zoom In Zoom Out Reset image size

During the second step, the element connectivity is not changed, solely the nodes are displaced by smoothing the mesh. Therefore, the well-known Laplacian mesh smoothing operation is applied [24]. In the algorithm, the domain vertex nodes are fixed. These nodes are not displaced during the entire algorithm, as it changes the domain geometry. The Laplacian smoothing step is applied first to nodes on a domain edge. The new positions of these nodes on the boundary are only influenced by the domain vertices and domain edge nodes. After this smoothing step, the domain edge nodes are temporarily fixed and the domain face nodes are displaced by the smoothing step. Finally, the interior nodes are displaced. In the Laplacian mesh smoothing operation, the new nodal position ${\overrightarrow{y}}_{p}^{\enspace \star }$ depends on old nodal positions of a set of neighboring nodes ${\overrightarrow{y}}_{l}$ according to

$$\begin{equation}{\overrightarrow{y}}_{p}^{\enspace \star }=\frac{1}{{N}_{l}}\sum\limits _{l\in {L}_{p}}\enspace {\overrightarrow{y}}_{l},\end{equation} \tag{ 12 }$$

in which N_l is the number of local nodes in the set L_p . Furthermore, a distinction is thus made between the handling of a given node p which is classified as domain vertex (B_p = 3), domain edge node (B_p = 2), domain face node (B_p = 1) or interior node (B_p = 0), as discussed before. Let A_p be all adjacent nodes of node p, then L_p is defined as

$$\begin{equation}{L}_{p}=\left\{n\in {A}_{p},{B}_{n}\geqslant {B}_{p}\right\}.\end{equation} \tag{ 13 }$$

By performing the smoothing operation in this manner, there are no special constraints needed for the displacement of nodes at the domain boundary if the domain shape is cuboid, which is the case for the example domain. However, in case of non-straight boundaries, a gradient descend algorithm is used to keep the nodes on the corresponding external boundaries. Therefore, the residual R_p of node p is defined as

$$\begin{equation}{R}_{p}=\sum\limits _{k\in {K}_{p}}\enspace {f\hspace{-2pt}}_{k}({\overrightarrow{y}}_{p}^{\enspace \star }),\end{equation} \tag{ 14 }$$

which vanishes when the node is located on the corresponding external boundary. The iterative update of nodal positions is then given by

$$\begin{equation}\delta {\overrightarrow{y}}_{p}^{\enspace \star }=-\frac{{R}_{p}}{(\overrightarrow{\nabla }{R}_{p}\cdot {\overrightarrow{e}}_{R})}{\overrightarrow{e}}_{R},\end{equation} \tag{ 15 }$$

with

$$\begin{equation}{\overrightarrow{e}}_{R}=\frac{\overrightarrow{\nabla }{R}_{p}}{{\Vert}\overrightarrow{\nabla }{R}_{p}{\Vert}}\end{equation} \tag{ 16 }$$

the unit vector in the direction of steepest ascend. This algorithm is iterated until the residual R_p drops below a given tolerance . Notice that in theory the residual R_p in equation (14) can also vanish if two boundary functions cancel each other. If this is the case for a specific application, the boundary functions ${f\hspace{-2pt}}_{k}(\overrightarrow{y}\enspace )$ can easily be adapted to represent the squared or unsigned distance to a domain boundary face, leading to ${f\hspace{-2pt}}_{k}(\overrightarrow{y}\enspace )\geqslant 0,\enspace \forall \enspace \overrightarrow{y}$. However, it should be noted that for these boundary functions, the unit vector in the direction of steepest ascend ${\overrightarrow{e}}_{R}$, see equation (16), becomes inaccurate at the domain boundary. This is the reason signed distance boundary functions are currently adopted.

Due to these local operations, the element quality of the pillow elements increases. The entire smoothing operation is iterated five times, which yields a good balance between obtaining a smooth grain conforming mesh and preserving the main grain shape features. In figure 4(b), the example case mesh is shown after the smoothing operations have been applied.

In figure 5(a), only the red and green grain of the example case are shown. Here, in the highlighted yellow region, an element with two edges on an internal junction between three grains (the red, green and blue grain, indicated by the blue line) is identified. Such elements are undesirable since smooth grain boundaries at these locations would introduce highly distorted or even invalid elements. Therefore, in the third step of the conversion to a grain conforming mesh, it is ensured that each element has maximum one edge on an internal junction between three or more grains. This is accomplished by another pillowing step (which is referred to as the surface pillowing step), in which a pillow layer is introduced around the set of elements that have one of their faces on a given internal grain boundary face. The rest of the pillowing operation is identical to the first step, which has been explained before. For the example case, this leads to the mesh as visualized in figure 4(c).

Zoom In Zoom Out Reset image size

The fourth and final step of the conversion to a grain conforming mesh is the same as the second step. The mesh is smoothened locally to increase the mesh quality of the new surface pillow elements and their neighbors. Figure 5(b) reveals that after the third and fourth step, each element has maximum one edge on a junction between three or more grains, allowing to obtain a fully grain conforming mesh, which, for the example case, is shown in figure 4(d).

For the example case, the mesh is now complete. It does not contain any invalid or highly distorted elements and the discretization can directly be used in microstructural simulations. However, for more complex cases (including the application cases described in section 4), it is necessary to increase the mesh quality or even eliminate inverted hexahedral elements. Therefore, the algorithm described in [25], which is specifically applicable to hexahedral meshes, is adopted. This algorithm deploys a nodal quality metric based on a normalized Jacobian. Given a hexahedral element e with node p and neighboring element nodes A, B and C (which are also a vertex of the same element e), as visualized in figure 6, let $\overrightarrow{a}$, $\overrightarrow{b}$ and $\overrightarrow{c}$ denote the edge vectors from node p to nodes A, B and C, respectively. Note that the ordering of these edge vectors is defined such that the volume $J=\overrightarrow{a}\cdot (\overrightarrow{b}\times \overrightarrow{c})$ spanned by these vectors is positive for a well-defined hexahedral element e. The quality ${Q}_{p}^{e}$ of node p with respect to element e is then given by

$$\begin{equation}{Q}_{p}^{e}=\hat{\overrightarrow{a}}\cdot (\hat{\overrightarrow{b}}\times \hat{\overrightarrow{c}}),\end{equation} \tag{ 17 }$$

in which

$$\begin{equation}\hat{\overrightarrow{a}}=\frac{\overrightarrow{a}}{{L}_{\mathrm{max}}},\qquad \hat{\overrightarrow{b}}=\frac{\overrightarrow{b}}{{L}_{\mathrm{max}}},\qquad \hat{\overrightarrow{c}}=\frac{\overrightarrow{c}}{{L}_{\mathrm{max}}},\end{equation} \tag{ 18 }$$

with L_max the maximum local element edge length, i.e.

$$\begin{equation}{L}_{\mathrm{max}}=\mathrm{max}\left\{{\Vert}\overrightarrow{a}{\Vert},{\Vert}\overrightarrow{b}{\Vert},{\Vert}\overrightarrow{c}{\Vert}\right\}.\end{equation} \tag{ 19 }$$

By defining the quality ${Q}_{p}^{e}$ in this manner, not only the orthogonality of element e matters, but also its aspect ratio. Furthermore, the quality Q_p of node p can now be defined as

$$\begin{equation}{Q}_{p}=\underset{e}{\mathrm{min}}\enspace {Q}_{p}^{e},\end{equation} \tag{ 20 }$$

representing the minimum quality ${Q}_{p}^{e}$ of all elements e that share node p. Finally, the mesh quality Q_mesh is defined as

$$\begin{equation}{Q}_{\text{mesh}}=\underset{p}{\mathrm{min}}\enspace {Q}_{p}.\end{equation} \tag{ 21 }$$

An invalid mesh can then be identified as Q_mesh ⩽ 0, indicating that the discretization contains at least one element with a negative volume. Furthermore, the mesh quality is considered poor if Q_mesh is smaller than a given minimum mesh quality Q_min.

Zoom In Zoom Out Reset image size

In order to improve the mesh quality Q_mesh, nodes are displaced. In general, multiple iterations are required to reach the minimum mesh quality Q_min. The condition to accept a new nodal position ${\overrightarrow{y}}_{p}^{\enspace \star }$ is therefore based on a positive iterative difference in quality ΔQ_p (and not on the absolute value of the nodal quality Q_p ). Furthermore, it is required that the new nodal position ${\overrightarrow{y}}_{p}^{\enspace \star }$ does not decrease the quality of any of its adjacent nodes A_p , i.e.

$$\begin{equation}{\Delta}{Q}_{p} > 0\quad \wedge \quad {\Delta}{Q}_{a}\geqslant 0\quad \forall \enspace a\in {A}_{p}.\end{equation} \tag{ 22 }$$

The nodes are sorted based on their quality Q_p and the node with the lowest quality is improved first. The algorithm uses two techniques to increase the nodal quality of node p, which are the same as described by Calvo and Idelsohn [25]. For the sake of completeness, they are briefly described below:

• (a)  
  Gradient driven.The update of the nodal position ${\Delta}{\overrightarrow{y}}_{p}$ is in the direction of the gradient of Q_p with a step size equal to a fraction ξ of the average local mesh size ${L}_{\text{avg}}=\frac{1}{3}({\Vert}\overrightarrow{a}{\Vert}+{\Vert}\overrightarrow{b}{\Vert}+{\Vert}\overrightarrow{c}{\Vert})$, i.e.

  $$\begin{equation}{\Delta}{\overrightarrow{y}}_{p}=\xi {L}_{\text{avg}}\frac{\overrightarrow{\nabla }{Q}_{p}}{{\Vert}\overrightarrow{\nabla }{Q}_{p}{\Vert}}.\end{equation} \tag{ 23 }$$

  The fraction ξ is iteratively decreased from a maximum of ξ = 1 to a minimum of $\xi =\frac{1}{128}$, each iteration halving the fraction. The first nodal displacement ${\Delta}{\overrightarrow{y}}_{p}$ to be accepted (i.e. condition (22) is fulfilled), is applied.
• (b)  
  Randomized.In case the first technique did not result in an accepted nodal displacement, the second technique is applied. Here a random vector on the unit sphere ${\overrightarrow{e}}_{\mathrm{s}}$ is adopted as the direction vector, and a random fraction 0 < ξ < 1 defines the step size, leading to

  $$\begin{equation}{\Delta}{\overrightarrow{y}}_{p}=\xi {L}_{\text{avg}}{\overrightarrow{e}}_{\mathrm{s}}.\end{equation} \tag{ 24 }$$

  The number of repetitions of taking a random direction vector ${\overrightarrow{e}}_{\mathrm{s}}$ and a random fraction ξ is linearly scaled with the quality Q_p , between 50 for Q_p ⩽ 0 and 0 for Q_p = 1. Here, also, the first nodal displacement ${\Delta}{\overrightarrow{y}}_{p}$ to fulfill condition (22), is applied.

In this manner, nodes in the interior of the domain (B_p = 0) are updated. Nodes on a domain boundary face (B_p = 1) require extra consideration. These are kept on the corresponding boundary by applying the same gradient descend algorithm as described in section 3.3. Naturally, the evaluation of the new nodal position ${\overrightarrow{y}}_{p}^{\enspace \star }$ is performed after convergence of the gradient descend iterations. Furthermore, it has been observed that it is not required for the nodes on a domain boundary edge (B_p = 2) to be moved, to obtain a valid mesh.

After the nodes are displaced individually, the new mesh quality Q_mesh is calculated. This entire procedure is performed iteratively until a given minimum mesh quality Q_min is reached. In order to speed up this mesh optimization algorithm, a node p for which Q_p > 2Q_min is skipped in the displacement update, since its update has a negligible effect on the new mesh quality Q_mesh. Nodes p with a quality Q_min < Q_p < 2Q_min are still updated using this procedure, since this update typically allows for more freedom of neighboring nodes with a lower quality.

The mesh optimization strategy as discussed is applied twice in the 3D anisotropic Voronoi algorithm, after both step 2 and step 4 of the conversion to a grain conforming geometry. It is thus used as an extra optimization step after the Laplacian mesh smoothing operation to ensure a good final mesh quality. As explained before, at the end of step 2 the element connectivity is not yet suitable for smooth internal junctions between three or more grains. Therefore, in the first mesh optimization step, nodes which are located on such edges are not taken into account for the determination of the mesh quality Q_mesh.

In micromechanical modeling applications, it is often desired to have a periodic discretization. This can also be achieved with the 3D anisotropic Voronoi algorithm using the general approach of generating periodic computational microstructures. However, a few important considerations should be highlighted.

First, the initial mesh is required to be periodic. In the algorithm, nodes from the initial mesh on opposite domain boundary faces, edges and vertices are identified and tied together. This also holds for the newly generated nodes during the pillowing operations. The ties are applied specifically after the Laplacian mesh smoothing operations and during the mesh quality improvement procedure. Here, it is also taken into account in the acceptance of a new nodal position.

Secondly, all grain seed points are copied 3 × 3 × 3 = 27 times in order to generate a periodic geometry in the center, but the initial mesh is not replicated to save computational costs. In the iterative acceptance of the assignment of elements to grains (section 3.2), the periodicity is taken into account in the definition of neighboring elements. Furthermore, in order for the pillowing operations to be periodic (and to introduce the new nodes periodically), the intersection line of a given grain boundary with a domain boundary face should be at the same location on opposite periodic boundary faces. Therefore, it is important that two periodic neighboring elements sharing a common element face on an opposite periodic domain boundary are assigned to the same grain. Equivalently, a set of elements sharing an edge (or vertex) on corresponding periodic domain edges (or corners), are assigned to the same grain. These considerations ensure that a periodic hexahedral mesh is created.

The first application case consists of an AISI 316LSi stainless steel microstructure which has been produced using wire + arc additive manufacturing (WAAM). This process is similar to multi-pass fusion welding. By stacking weld lines on top of each other, a large-scale metallic component (∼1 to $>$10 m) can be built up. Locally, an electric arc is used to melt the material during deposition. When this melt pool solidifies, a new bead of solid metal results. During solidification, crystals tend to grow in the direction of the local temperature gradient, which is strongly present during the manufacturing process. This promotes the growth of strongly columnar grains and, due to the melt pool shape, the columnar direction is spatially varying.

Electron backscatter diffraction (EBSD) data of a typical grain structure for austenitic stainless steels produced by the WAAM process obtained from [33] is depicted for two cross-sectional planes: figure 7(a) shows the deposition direction—building direction plane and figure 7(b) displays the transverse direction—building direction plane, which have been generated using MTEX [34]. Here, the deposition direction is the direction in which the electric arc moves, the building direction is the vertical stacking direction of layers and the transverse direction is the direction of neighboring beads in a single layer. In order to generate a 3D statistically representative and periodic computational microstructure, the 3D anisotropic Voronoi algorithm is applied. The domain of interest and corresponding boundary functions are acquired from the repeated stacking of weld lines. An experimentally obtained average 2D periodic fusion zone shape [35], which is non-convex and contains curved boundaries, is extruded in the third dimension to account for about four grains in the deposition direction, which is based on the average grain thickness obtained from figure 7(a). The resulting domain edges are visualized in figure 8(a) by black lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The grain seed points used as input for the algorithm are obtained from the EBSD data displayed in figure 7(b). For this purpose, first, grains are identified from the EBSD data by a misorientation threshold of 5°. Ellipses are fitted to the grains by equalizing their second moments [36]. The ellipse centroids obtained from the four regions in figure 7(b) are aligned according to their corresponding experimental fusion zone boundaries estimated from the EBSD data and stacked upon each other in the third (deposition) direction. The ellipses are symmetric around their long axes to form a spheroid. This assumes that the short axis of the grain in the transverse direction—building direction cross-section is the same as the thickness of the grain in the deposition direction, which is supported by the EBSD data shown in figure 7(a). Finally, a random normal distribution with a standard deviation of 10% of the layer height is applied to the grain seed point coordinates in the third (deposition) direction. This then results in the final grain seed point input for the algorithm as displayed in figure 8(a). Here, each spheroid represents the velocity field of an individual grain seed point. Their color represents the grain orientation with respect to the building direction, which is obtained from the corresponding grain in the EBSD data.

As mentioned before, the domain of interest and corresponding boundary functions are acquired from the repeated stacking of weld lines. The experimentally obtained average 2D periodic fusion zone shape [35] is discretized into quad elements and subsequently extruded to form 10 layers of hexahedral elements in the out-of-plane (deposition) direction. Although this is a relatively coarse initial discretization in the deposition direction, in the final mesh, each grain consists of at least 5 elements in their smallest dimension, due to the pillow element insertion step. The initial mesh consists of 54 100 hexahedral elements and 61 171 nodes as shown in figure 8(b). Here, the state of the algorithm after the first step (assignment of elements to grains) is illustrated by the element coloring.

The final computational geometry and mesh, generated by the 3D anisotropic Voronoi algorithm, are depicted in figures 8(c) and (d), respectively. The former shows a 3D computational microstructure consisting of smooth grain boundaries and strongly columnar grains with their long axes corresponding to the experimental data shown in figure 7(b). The latter consists of a periodic hexahedral and grain conforming mesh (381 616 elements and 415 839 nodes) that can be used directly in e.g. finite element crystal plasticity simulations.

In order to analyze the representativeness of the computational microstructure with respect to the experimental data, the same cross-sections on which the EBSD data (figure 7) is gathered, are obtained from the final computational microstructure (figure 8(c)). These cross-sections are shown in figure 9, in which the color represents the grain orientation in the building (vertical) direction. During the following analysis, the periodicity of the microstructure is taken into account. Therefore, grains that continue over the domain boundary are treated as a single grain. This leads to the grain shapes illustrated in figure 9 by black grain boundaries. In these figures, the light-shaded grains in the background serve to emphasize the periodicity of the microstructure and these are not taken into account in the analysis.

Zoom In Zoom Out Reset image size

A comparison is made in terms of the grain size distribution and columnar grain direction and aspect ratio, for both cross-sectional planes, to evaluate the correspondence between the EBSD data and the resulting computational microstructure. First, focus is put on the deposition direction—building direction plane. Figure 10(a) shows the total percentage of area (cumulative grain size) that is covered by grains having a certain grain size or smaller in this cross-section. Depending on the choice of misorientation threshold used in the definition of a sub-grain ($\sim 0.5^{\circ}$) or a grain ($\sim 5^{\circ}$) in the EBSD data, a range of grain size distributions can be obtained. This is illustrated by the gray shaded area representing misorientations between 0.5°–5° as indicated. Obviously, the resulting computational microstructure adequately represents the experimental grain size distribution. The grain columnar direction and aspect ratio are displayed in figure 11(a), where each dot represents a single grain. Both quantities are defined with respect to a fitted ellipse: the former defines the orientation of the major axis and the latter specifies the ratio between the major and minor axis length. Also here, the resulting computational microstructure represents the EBSD data appropriately. Notice that the large grain size variations in the experimental data, see figure 10(a), are not obtained with the computational microstructure due to the limited thickness of the domain and the small number of grains in the deposition direction. Also, due to the relatively large EBSD scan with many grains in this cross-section, a wider spread of experimental data is observed in both figures 10(a) and 11(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figures 10(b) and 11(b), the same comparison is shown for the transverse direction—building direction plane. Note that here, the data of the four cross-sections shown in figures 7(b) and 9(b) are combined. Again, the computational microstructure generated by the 3D anisotropic Voronoi algorithm adequately represents the experimental EBSD data. Some minor deviations can be observed in figure 10(b) in the representation of small grains. This can be attributed to the initial mesh size. Some grain seed points with a very small volume have not been assigned to an element in the initial mesh. Upon refinement, more small grains can be taken into account. Such refinement would also improve the resulting grain aspect ratios. The insertion of a total of four pillow elements over the minor axis of an elongated grain has a much larger effect on the minor axis length than the effect of insertion of the same four pillow elements over the major axis of the same grain on the major axis length. Therefore, due to the pillow element insertion step of the algorithm, the aspect ratio of a largely elongated grain decreases slightly. A finer initial discretization reduces this effect and improves the representativeness of the computational microstructure even further.

To show the wide range of applicability of the algorithm, the second application case consists of a tension loaded magnesium AZ31B alloy with a twinned microstructure. Although twinning is not the mechanism that is mimicked in the 3D anisotropic Voronoi algorithm, the resulting experimental microstructure consists of elongated grains with columnar directions that vary spatially, which can still be represented by the algorithm. 3D EBSD data of a magnesium alloy microstructure is obtained from [28] and a subsection of 18.4 × 20.8 × 20 μm is shown in figure 12(a).

Zoom In Zoom Out Reset image size

The input of the 3D anisotropic Voronoi algorithm is obtained by fitting ellipsoids to grains by again equalizing their second moments. In this case, some grains have been fitted by multiple ellipsoids through first splitting each grain into various bulky parts. However, note that the same orientation is assigned to each grain part, so that it becomes a single grain in the final geometry. Each of these regions is represented by a single ellipsoid. Its center is subsequently used as seed point coordinate and the ellipsoid shape represents the corresponding velocity growth field. The resulting input of the algorithm is illustrated in figure 12(b).

A voxel mesh of 16 × 18 × 17 = 4896 hexahedral elements is adopted as initial mesh, which is depicted in figure 12(c), together with the algorithm's assignment of elements to grains. Periodic boundary conditions are not applied, since the purpose of this application case is to represent a specific part of the magnesium alloy microstructure. The final computational geometry and corresponding hexahedral mesh are shown in figures 12(d) and (e), respectively. Upon comparison of the experimental 3D EBSD data (figure 12(a)) with the final computational grain geometry (figure 12(d)), a convincing resemblance is observed. The final hexahedral mesh is fully grain conforming and consists of 27 964 elements and 31 626 nodes. Note that using this discretization (94 878 degrees of freedom) in microstructural simulations is much more computationally efficient compared to using the experimental 3D EBSD data as direct voxel mesh (381 123 degrees of freedom). On top of that, using the experimental data directly, results in the presence of non-grain-conforming boundaries, possibly introducing spurious stress concentrations at these staggered interfaces [21], which are not present in the final mesh of the 3D anisotropic Voronoi algorithm.

In order to analyze the correspondence in more detail, a comparison is made between the grain seed point input data of the algorithm and the characteristics of the output grain morphology in terms of grain size, aspect ratio and orientation. In figure 13(a), the input and resulting grain size distribution are compared. A convincing correspondence is observed. In figure 13(b), the aspect ratio (maximum over minimum principal axis length) of the grains in the input and result of the algorithm are compared. Also here, as satisfactory similarity is established. Similar to the WAAM computational microstructure, the anisotropic 3D Voronoi algorithm makes largely elongated grains slightly more spherical, due to the pillow element insertion step.

Zoom In Zoom Out Reset image size

In the pole figures presented in figures 14(a) and (b), each dot represents a single normalized principal axis direction vector of the input velocity field and the resulting grain fitted ellipsoid of the computational microstructure, respectively. The color of each dot resembles the ratio of the corresponding principal axis length to the minimum principal axis length for each ellipsoid. Focusing on the positions of the dots, it is obvious that the general characteristics of the principal axis directions are reproduced properly by the 3D anisotropic Voronoi algorithm. The colors of the dots also reveal that the procedure slightly decreases the aspect ratio of largely elongated grains. If needed, this effect can be reduced by employing a finer initial mesh to improve the representativeness of the computational microstructure even further.

Zoom In Zoom Out Reset image size