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

At the microstructural scale, Voronoi tessellations are commonly used to represent a polycrystalline morphology. However, due to spherical growth of nuclei, an anisotropic tessellation with spatially varying elongated grain directions, which is present in many applications, cannot be obtained. In this work, a novel 3D anisotropic Voronoi algorithm is presented, together with its implementation and two application cases. The proposed algorithm takes into account preferred grain growth directions, aspect ratios and sizes in the definition of an ellipsoidal growth velocity field defined per grain. For applications in which a predetermined mesh is used, e.g. voxel-mesh based simulations, the grains are extracted in a straight-forward manner. In cases where a fully grain conforming discretization is desired, e.g. finite element simulations, a hexahedral mesh generator is incorporated to arrive at a discretization which can be directly used in microstructuralmodeling simulations. Two application cases are studied (a wire + arc additively manufactured and a magnesium alloy microstructure) in which the algorithm’s capability for curved, non-convex, periodic domains is shown. Furthermore, the resulting grain morphology is compared to experimental data in terms of grain size, grain aspect ratio and grain columnar direction distribution. In both cases, the algorithm adequately produces a representative volume

At the microstructural scale, Voronoi tessellations are commonly used to represent a polycrystalline morphology.However, due to spherical growth of nuclei, an anisotropic tessellation with spatially varying elongated grain directions, which is present in many applications, cannot be obtained.In this work, a novel 3D anisotropic Voronoi algorithm is presented, together with its implementation and two application cases.The proposed algorithm takes into account preferred grain growth directions, aspect ratios and sizes in the definition of an ellipsoidal growth velocity field defined per grain.For applications in which a predetermined mesh is used, e.g.voxel-mesh based simulations, the grains are extracted in a straight-forward manner.In cases where a fully grain conforming discretization is desired, e.g.finite element simulations, a hexahedral mesh generator is incorporated to arrive at a discretization which can be directly used in microstructural modeling simulations.Two application cases are studied (a wire + arc additively manufactured and a magnesium alloy microstructure) in which the algorithm's capability for curved, non-convex, periodic domains is shown.Furthermore, the resulting grain morphology is compared to experimental data in terms of grain size, grain aspect ratio and grain columnar direction distribution.In both cases, the algorithm adequately produces a representative volume element with convincing representativeness of the experimental data.The 3D anisotropic Voronoi algorithm is highly versatile in a wide range of application

Introduction
In many engineering applications, the macroscopic properties of materials are predicted based on detailed modeling of the underlying microstructure.Being able to generate specific geometrical features of the microstructure is thereby an important ingredient.Voronoi tessellations are frequently used to represent grain shapes of (metal) polycrystalline materials at the microstructural scale [1][2][3].Such tessellations are generated starting from a set of grain seed points, which represent the grain nuclei.In the Voronoi tessellation, a spherical growth from these seed points outward, is assumed [4].The interface where two growth fronts meet, represents a grain boundary.As a result, a specific grain is defined by the volume that is closer to the corresponding grain seed point than to all other grain seed points.
For homogeneously distributed grain seed points, an isotropic grain morphology is obtained, due to the assumption of spherical growth.However, for many applications, this does not properly represent the actual microstructure.Therefore, many extensions to the standard Voronoi algorithm have been used to improve the microstructural description.This includes e.g.hardsphere Voronoi [5,6], Laguerre-Voronoi [7,8] or the Johnson-Mell tessellation [9,10].Using such algorithms, grain size distributions, curved and concave-shaped grains can be incorporated.However, all these extensions still rely on an isotropic (spherical) growth speed.This is not representative for applications in which columnar grains with (spatially varying) columnar directions are present.This typically applies to microstructures obtained by additive manufacturing [11,12], casting [13], welding [14], etc, i.e. methods which have a 'directional' driving force.For the modeling of such microstructures, the standard Voronoi approach would not reproduce the correct grain morphology since a spherical growth velocity field is not appropriate.
In previous work [15], an anisotropic Voronoi algorithm was developed which generates two-dimensional computational microstructures based on an elliptical grain growth velocity field to incorporate spatially varying grain growth directions.The output of this algorithm consists of (convex) grain vertices and connectivity.Finite element discretization of the geometry is applied afterwards.In this work, the algorithm is extended to three dimensions by considering an ellipsoidal grain growth velocity field and by including the mesh generation step.This enables the generation of concave shaped grains, increasing the representativeness of the generated microstructure.
Various standard element types are generally available in mesh generation algorithms and in commercial finite element packages.Although in 3D, tetrahedral meshes, based on Delaunay triangulation [16], are generally easier to generate compared to hexahedral meshes, tetrahedral meshes also bring serious disadvantages [17].Especially, tetrahedral meshes tend to show volumetric locking (artificial stiffening) when used in finite element simulations of (nearly) incompressible materials, such as based on crystal plasticity, requiring the use of e.g.nonstandard tetrahedral elements [18].Therefore, in this work, specific attention is given to the generation of an all-hexahedral based mesh, which can be directly used for microstructural simulations.
In general, a distinction can be made between geometry-first and mesh-first finite element discretization procedures [19].In the former, the surface of the domain is used as a starting point for the mesh to develop inwards.The latter uses a simple (or even regular) mesh that covers the entire domain of interest.Then, this initial mesh is modified locally to capture specific features [20].The 3D anisotropic Voronoi algorithm generates grains and a final grain-conforming mesh is essential to make sure no spurious stress concentrations arise at staggered, non-grain-conforming boundaries [21].Hence, a mesh-first generation procedure is adopted in the current work.
This paper is organized as follows.In section 2, the mathematical definition of the ellipsoidal grain growth velocity field is specified.Then, in section 3, the numerical implementation of the algorithm is detailed.Furthermore, the 3D anisotropic Voronoi algorithm is applied to two cases in section 4, along with an analysis of the generated microstructure correspondence to input data.Finally, conclusions are drawn in section 5.

Mathematical background
The definition of a grain C i is based on the time t i required to grow from the grain seed point x i to an arbitrary point in the domain of interest y, with respect to a Cartesian coordinate system { e 1 , e 2 , e 3 }, as in which i, j = 1, . . ., N, with N the number of grain seed points.The time t i is obtained from in which d( a, b) ≡ a − b is the Euclidean distance function and v i ( y ) is the (scalar) growth speed in the direction y − x i .
In order to take into account the preferred growth directions of the grains, an ellipsoidal velocity field v i ( y) is now introduced to describe the anisotropic growth speed of each grain.Given grain seed point i, an arbitrary ellipsoidal velocity field can be defined by the volume of the ellipsoid V i and two preferred growth vectors, p i and q i , which are perpendicular to each other (i.e.p i • q i = 0).The directions and lengths of these vectors define the principal axes and corresponding aspect ratios of the ellipsoid, respectively, as defined in appendix A and as depicted in figure 1.The resulting description of the directional growth speed v i can be expressed in terms of the directional unit vector e r = e r ( y ) = y− x i y− x i as in which the scaling factor s i is defined by and r3 is the dimensionless radius of the 3D ellipsoid defined as r3 r2 ( e r ) 2 + with The derivation of these equations is given in appendix A. In these equations, the principal directions e p = p i p i and e q = q i q i are used.Furthermore, r2 represents the dimensionless radius of the 2D elliptical cross-section of the 3D ellipsoid in the plane through the center and orthogonal to q i , which is highlighted in figure 1.Consequently, e r is the normalized projection of e r onto this plane, i.e. e r = e r − e r • e q e q e r − e r • e q e q .Note that in case the anisotropic growth field is spheroidal (i.e. an ellipsoid with two axes of the same length), the ellipsoid is rotationally symmetric around its single preferred growth direction p i and the description of the velocity field v i reduces to

Input
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 1 ( y), . . ., f N b ( y ) and their gradients ∇ f 1 ( y ), . . ., ∇ f N b ( y ) whose definition will be provided later, (iii) the seed point coordinates x 1 , . . ., x N , (iv) the corresponding preferred growth direction vectors p 1 , . . ., p N (including optional second preferred growth vectors q 1 , . . ., q N ) and (v) the corresponding volumes of the ellipsoids V 1 , . . ., 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 k ( y) should be provided.Each of these functions f k ( y ) are defined to be zero at the kth face of the domain, Γ k , reading 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 y p , let K p contain the set of domain boundary face numbers k on which this node is located, i.e.
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 ( x i , p i , 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 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 × 1 × 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 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.

Assignment of elements to grains
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 y = y e is the center of element e.Based on these growth times t 1 ( y e ), . . ., t N ( 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.
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 (j) after finishing.In (i) only the elements corresponding to the blue grain are shown.Note, the light-shaded color of the elements corresponds to the color of the preferred grain seed point assignment in figure 2. A dark-shaded element corresponds to the actually assigned grain seed point.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 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 nonmanifold conditions [22].For each node p of element e, the Euler number E e p,i 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 e p,¬i (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 e p,i = 2 or E e p,¬i = 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.

Converting to a grain-conforming mesh
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 nongrain-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 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 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.
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 y p depends on old nodal positions of a set of neighboring nodes y l according to 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 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 which vanishes when the node is located on the corresponding external boundary.The iterative update of nodal positions is then given by with 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 k ( y ) can easily be adapted to represent the squared The blue line indicates the internal junction of the red, green and blue grain (of which the latter is not shown).Furthermore, the yellow circle highlights the influence of step 3 and 4 on the grain conformity of the mesh.
or unsigned distance to a domain boundary face, leading to f k ( y ) 0, ∀ y.However, it should be noted that for these boundary functions, the unit vector in the direction of steepest ascend 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).
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).

Increasing mesh quality
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 a, b and 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 = a • ( b × c) spanned by these vectors is positive for a well-defined hexahedral element e.The quality Q e p of node p with respect to element e is then given by in which with L max the maximum local element edge length, i.e.
By defining the quality Q e p 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 representing the minimum quality Q e p of all elements e that share node p. Finally, the mesh quality Q mesh is defined as 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 .
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 y p 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 y p does not decrease the quality of any of its adjacent nodes A p , i.e.
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 Δ 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 The fraction ξ is iteratively decreased from a maximum of ξ = 1 to a minimum of ξ = 1 128 , each iteration halving the fraction.The first nodal displacement Δ 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 e s is adopted as the direction vector, and a random fraction 0 < ξ < 1 defines the step size, leading to The number of repetitions of taking a random direction vector e 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 Δ 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 y p 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 .

Generation of a periodic 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.

Application
The algorithm to generate 3D anisotropic Voronoi microstructures along with their finite element discretization can be applied to generate computational microstructures in a variety of fields, and is specifically suitable to polycrystalline microstructures in which grains with spatially varying elongated shape directions are a typical feature.This includes additive manufacturing [12,26,27], twinning [28,29], casting [13], welding [30,31], (dynamic) recrystallization [32], etc.In this section, the algorithm is applied to two application cases.In the first case, the algorithm's capability of generating a statistically representative volume element (RVE) is demonstrated for a unit cell with non-straight periodic boundaries.In the second application case, 3D EBSD data is directly used as input for the algorithm to generate a mesh that can be used to model a specific part of the microstructure.In both cases, grains with spatially varying elongated directions are characteristic features of the selected microstructures.

Wire + arc additive manufacturing
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.
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.
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 (∼0.5 • ) or a grain (∼5 • ) 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).
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.Comparison of the grain size distribution between the experimental EBSD data and the resulting computational microstructure for the WAAM example in two cross-sectional planes: (a) deposition direction-building direction and (b) transverse direction-building direction.The gray shaded area illustrates the range of experimental grain size distributions that can be obtained for a grain misorientation threshold taken between 0.5 • -5 • as indicated.Comparison between the experimental EBSD data for two grain misorientation threshold definitions (as indicated) and resulting computational geometry for the WAAM microstructure in two cross-sectional planes: (a) deposition direction-building direction and (b) transverse direction-building direction.Each dot represents a single grain.The angular coordinate corresponds to its columnar direction and the radial coordinate corresponds to its aspect ratio.

Magnesium alloy
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).
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.Equal area pole figures in which each dot represents a single normalized principal axis direction vector of (a) the input velocity field and (b) the grain fitted ellipsoids of the computational geometry.The color of each dot resembles the ratio of the corresponding principal axis length to the minimum principal axis length for each ellipsoid.This ratio trivially equals 1 for the principal axis vector with minimum length, which is not relevant.These vectors are therefore depicted by an open circle.
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.
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.

Conclusions
In this work, a previously developed two-dimensional anisotropic Voronoi algorithm is extended to 3D, along with the generation of a hexahedral mesh and mesh optimization.The algorithm takes into account the preferred grain growth directions, aspect ratios and sizes in the definition of an ellipsoidal growth velocity field, which is used to assign elements from a simple initial hexahedral mesh with corresponding boundary definitions to specific grain seed points.For specific applications in which a predetermined mesh is used, e.g. using an FFT solver, the assignment of elements to grain seed points can be considered the output of the algorithm.For many other applications, e.g.finite element crystal plasticity simulations, element pillowing, smoothing and mesh optimization operations are applied afterwards to arrive at a fully grain conforming and hexahedral mesh, which can be directly used in microstructural modeling simulations.The required input for the algorithm can be obtained from e.g.thermal (history) simulations, microstructural statistics or (3D) experimental observations.
The algorithm is useful in a wide variety of applications, and is specifically interesting for the generation of computational microstructures that include grains with a spatially varying elongated geometry.Two application cases have been discussed, a wire + arc additively manufactured microstructure and a magnesium alloy with a twinned microstructure.It is demonstrated that the algorithm can be used in curved, non-convex domains and that periodicity of the grains and the mesh can be enforced.The computational microstructure is compared to experimental data in terms of grain size, grain aspect ratio and grain columnar direction distribution.In both cases, the algorithm adequately reproduces an RVE with a convincing resemblance to the experimental data.
From this equation, the directional growth velocity v( e r ) can be found as v( e r ) = r = sr 3 , (37) with normalized radius r3 defined as which is more useful than equation (33) in the numerical implementation.In this equation, e r = e r − e r • e q e q e r − e r • e q e q (40) is the normalized projection of e r onto the 2D elliptical cross-section of the 3D ellipsoid in the plane through the center and orthogonal to q (see figure 1).

Figure 1 .
Figure 1.Schematic representation of the 3D ellipsoidal growth velocity field.The highlighted blue cross-section represents the plane in which e r is defined.

Figure 2 .
Figure 2. Illustration of the example case input.

Figure 3 .
Figure 3. Stages of the acceptance algorithm for the assignment of elements to grains: (a) before the start; (b)-(e) after the 1st, 2nd, 3rd and 4th iteration, respectively; (f) after the 18th iteration; (g) before the 19th iteration; (h) and (i) after the 22nd iteration;(j) after finishing.In (i) only the elements corresponding to the blue grain are shown.Note, the light-shaded color of the elements corresponds to the color of the preferred grain seed point assignment in figure2.A dark-shaded element corresponds to the actually assigned grain seed point.

Figure 4 .
Figure 4. Steps involved in the conversion to a grain conforming mesh.The state of the example case is visualized after (a) step 1: volume element insertion, (b) step 2: Laplacian mesh smoothing, (c) step 3: surface element insertion and (d) step 4: Laplacian mesh smoothing.

Figure 5 .
Figure 5. Section of the example case after (a) step 2 and (b) step 4 of the conversion to a grain conforming mesh.All elements belonging to the red and green grain are depicted.The blue line indicates the internal junction of the red, green and blue grain (of which the latter is not shown).Furthermore, the yellow circle highlights the influence of step 3 and 4 on the grain conformity of the mesh.

Figure 6 .
Figure 6.Hexahedral element with node p and corresponding neighbors A, B and C. Note that the ordering of nodes A, B and C with respect to node p matters.Adapted from [25].

Figure 7 .
Figure 7. EBSD data of the WAAM microstructure in two cross-sectional planes: (a) deposition direction-building direction (6.5 × 6.5 mm) and (b) transverse direction-building direction [33].In (b) four different scans of the same cross-section are shown, which all have a size of 4.0 × 2.4 mm, which is approximately one periodic unit.Colors indicate crystallographic directions in the building direction.

Figure 8 .
Figure 8. Four stages of the 3D anisotropic Voronoi algorithm during the computational WAAM microstructure generation: (a) input and domain of interest, (b) initial mesh and assignment of elements to grains, (c) final geometry and (d) final hexahedral mesh.The domain of interest has a periodic size of 3.85 × 2.34 × 0.50 mm.Colors correspond to grain orientations with respect to the building direction.

Figure 9 .
Figure 9. Periodic cross-sections of the computational WAAM microstructure (figure 8(c)) corresponding to the same planes as used for the EBSD data (figure 7): (a) deposition direction-building direction and (b) transverse direction-building direction.The cross-section in (a) is taken at the location indicated by the vertical dashdot line in (b).The solid black lines indicate grain boundaries of the grains used in the analysis and the light-shaded background grains illustrate the microstructural periodicity.This is also emphasized by the dashed line, indicating the periodic domain boundary.Colors indicate crystallographic directions in the building direction.

Figure 10 .
Figure10.Comparison of the grain size distribution between the experimental EBSD data and the resulting computational microstructure for the WAAM example in two cross-sectional planes: (a) deposition direction-building direction and (b) transverse direction-building direction.The gray shaded area illustrates the range of experimental grain size distributions that can be obtained for a grain misorientation threshold taken between 0.5 • -5 • as indicated.

Figure 11 .
Figure 11.Comparison between the experimental EBSD data for two grain misorientation threshold definitions (as indicated) and resulting computational geometry for the WAAM microstructure in two cross-sectional planes: (a) deposition direction-building direction and (b) transverse direction-building direction.Each dot represents a single grain.The angular coordinate corresponds to its columnar direction and the radial coordinate corresponds to its aspect ratio.

Figure 12 .
Figure 12.(a) Experimental 3D EBSD data of the magnesium alloy [28] and four stages of the 3D anisotropic Voronoi algorithm during the computational microstructure generation: (b) input and domain of interest, (c) initial mesh and assignment of elements to grains, (d) final geometry and (e) final hexahedral mesh.The top and bottom rows show the front and back of the domain, respectively.

Figure 13 .
Figure 13.Comparison of the (a) grain size distribution and (b) aspect ratio (maximum over minimum axis length) distribution between the input and resulting computational geometry generated by the 3D anisotropic Voronoi algorithm.

Figure 14 .
Figure14.Equal area pole figures in which each dot represents a single normalized principal axis direction vector of (a) the input velocity field and (b) the grain fitted ellipsoids of the computational geometry.The color of each dot resembles the ratio of the corresponding principal axis length to the minimum principal axis length for each ellipsoid.This ratio trivially equals 1 for the principal axis vector with minimum length, which is not relevant.These vectors are therefore depicted by an open circle.
manner, it can be shown that r2 can be written asr2 ( e r ) = 1 + 1 p i 2 − 1 e p • e r