Modeling
Geometry
Grains in polycrystalline microstructures can be approximated by ellipsoids. To generate synthetic microstructures, packing the particles (ellipsoids) which follow a particular size distribution into a pre-defined domain becomes the objective. The general framework employed in this regard is the collision detection and response system for particles under random motion in the box. Each particle \(i\) in the domain is defined as an ellipsoid in the three-dimensional Euclidean space \(\mathbb{R}^3\), and random position and velocity vectors \(\mathbf{r}^i\), \(\mathbf{v}^i\) are assigned to it. During their motion, the particles interact with each other and with the simulation box. The interaction between particles can be modeled by breaking it down into stages of collision detection and response, and the interaction between the particles and the simulation box can be modeled by evaluating if the particle crosses the boundaries of the box. If periodicity is enabled periodic images on the opposite boundaries of the box are created, otherwise the particle position and velocity vectors have to be updated to mimic the bouncing back effect.
Two layered collision detection
For \(n\) ellipsoids in the domain, the order of computational complexity would be \(O(n^2)\), since each ellipsoid in the domain is checked for collision with every other ellipsoid. A two-layered collision detection scheme is implemented to overcome this limitation. The outer layer uses an Octree data structure to segregate the ellipsoids into sub-domains; this process is done recursively until there are only a few ellipsoids left in each sub-domain. The inner layer consists of a bounding spheres hierarchy, wherein ellipsoids of each Octree sub-domain are tested for collision only when their corresponding bounding spheres overlap. This effectively reduces the number of collision checks and thus the order of computational complexity to \(O(nlog(n))\). The general framework of collision detection response systems with two-layer spatial-partitioning data structures has the following features:
Recursively decompose the given domain into sub-domains based on the Octree data structure.
Perform collision tests between bounding spheres of ellipsoids belonging to the same sub-domain.
Test for ellipsoid overlap condition only if the bounding spheres overlap.
Update the position and the velocity vectors \(\mathbf{r}\) and \(\mathbf{v}\) based on the collision response.
Test for collision with the simulation domain and create periodic images on the opposite boundaries or mimick the bouncing against the wall effect.
Layer 1: Octree data structure
To ensure efficient collision checks an Octree data structure is initialized on the simulation box. With pre-defined limits for Octree sub-division and particle assignment, the Octree trunk gets divided into sub-branches recursively. Thus, by only performing collision checks between particles belonging to a particular sub-branch, the overall simulation time is reduced.
Figure: Simulation domain and its successive sub-branches on two levels, where particles are represented by red filled circles (left). Octree data structure depicting three levels of sub-divisions of the tree trunk (right).
Layer 2: Bounding sphere hierarchy
Figure: Upper layer consists of the Octree sub-branch with particles (left), and lower layer is defined by bounding spheres for particles \(i, j\) with radii \(a^i, a^j\) respectively (right).
Overlap detection
The actual overlap of two static ellipsoids is determined by the algebraic separation condition developed by Wang (2001). Consider two ellipsoids \(\mathcal{A}: \mathbf{X}^T \mathbf{A} \mathbf{X} = 0\) and \(\mathcal{B}: \mathbf{X}^T \mathbf{B} \mathbf{X} = 0\) in \(\mathbb{R}^3\), where \(\mathbf{X} = [x, y, z, 1]^T\), the characteristic equation is given as,
Wang et al. have established that the equation has at least two negative roots, and depending on the nature of the remaining two roots, the separation conditions are given as,
\(\mathbf{A}\) and \(\mathbf{B}\) are separated if \(f(\lambda) = 0\) has two distinct positive roots.
\(\mathbf{A}\) and \(\mathbf{B}\) touch externally if \(f(\lambda) = 0\) has a positive double root.
\(\mathbf{A}\) and \(\mathbf{B}\) overlap for all other cases.
Particle (Ellipsoid) packing
The user-defined simulation box size and the ellipsoid size distribution are used for creating the simulation box and ellipsoids. The simulation begins by randomly placing ellipsoids of null volume inside the box, and each ellipsoid is given a random velocity vector for movement. As the simulation proceeds, the ellipsoids grow in size along their axes and also collide with one another updating their position and velocities. The simulation terminates once all the ellipsoids have reached their defined volumes; the process is depicted pictorially in the figure below.
Figure: Ellipsoid packing simulation with partice interactions at three different timesteps.
Note
Since the application is microstructure generation, where all grains have a predefined tilt angle, the angular velocity vector \((\mathbf{w})\) is not considered for the ellipsoids and thus their orientations are constrained.
Texture
In this section a brief summary of the Orientation Distribution Function (ODF) reconstruction is presented. A detailed description of this can be found in Biswas (2020) as a \(L_1\) minimization scheme. Furthermore, an orientation assignment algorithm is presented which takes the grain boundary texture into consideration during the assignment process.
ODF reconstruction with discrete orientations
Crystallographic texture can be represented in the form of a continuous functions called ODF i.e., \(f:SO(3) \rightarrow \mathbb{R}\). With the availability of Electron Back Scatter Diffraction (EBSD) equipment, the polycrystalline materials are easily characterized in the form of measured crystallographic orientations \(g_i, \ i=1,...,N\).
The measured orientations \(g_i\) are used to estimate the ODF by a kernel density estimation Hielscher (2013). A bell shaped kernel function \(\psi: \ [0,\pi] \rightarrow \mathbb{R}\) is placed at each \(g_i\), which when combined together estimates the ODF as
where \(\omega(g_i , g)\) is the disorientation angle between the orientations \(g_i\) and \(g\). It is vital to note here that the estimated ODF \(f\) heavily depends on the choice of the kernel function \(\psi\). Keeping this in mind, the de la Vall'ee Poussin kernel \(\psi_\kappa\) is used for ODf reconstruction within Kanapy. Please refer to Schaeben (1997) for a detailed description of the de la Vall'ee Poussin kernel function as well as its advantages.
Within the numerical modeling framework for polycrystalline materials, the micromechanical modeling requires a reduced number \(N^\prime\) of discrete orientations \({g^\prime}_i\) to be assigned to grains in the RVE. The ODF of the reduced number of orientations is given as,
Since \(N^\prime \ll N\), the the kernel shape parameter \(\kappa\) must be optimized (\(\kappa^\prime\)) such that the ODF estimated from \({g^\prime}_i\) should be close to the input ODF \(f\). To quantify the difference between them an \(L_1\) error can be defined on the fixed \(SO(3)\) grid as
Orientation assignment process
In addition to the crystallographic texture, polycrystalline material also have grain boundary texture which is represented in the form of the misorientation distribution function (MDF). Similar to ODF, MDF can be estimated on \(\Delta \phi:SO(3) \rightarrow \mathbb{R}\) due to the disorientation (\(\Delta g_i, \ i=1,...,N_g\)) at grain boundary segments. These can be used to assign the discrete orientations obtained after ODF reconstruction to the grains in the RVE generated by Kanapy’s geometry module. Both the orientations and the disorientations play different roles in the mechanical behavior of the material, a detailed discussion of which can be found in Biswas (2019).
Other than the crystallographic orientation of the grains, the dimension of the grain boundary is also a key aspect. Therefore, to incorporate the effect of the grain boundary dimension, \(\Delta g\) at each segment is weighted as per the segment dimension, as suggested in Kocks (2000). Since the EBSD data is a 2D image this weighting factor is estimated as \(w_i = s_i/S\), where \(s_i\) is corresponding segment length and \(S = \sum_i^{N_g} s_i\).
The disorientation \(\Delta g\) can be represented in axis-angle notation. And within Kanapy’s algorithms we focus on the angle (\(\omega\)) part of \(\Delta g\), commonly referred to as the disorientation angle. To imitate the statistical distribution from experiments, a Monte-Carlo scheme is suggested in Miodownik (1999). This is used here in a \(L_1\) minimization framework similar to the ODF reconstruction discussed earlier.
The assignment algorithm begins by randomly assigning orientations obtained from ODF reconstruction to the grains in the RVE. The disorientation angle distribution is estimated in the present configuration including the weights (due to the corresponding grain boundaries) and the \(L_1\) error is estimated. The orientations are then exchanged between two grains modifying the configuration. The \(L_1\) error is estimated for the modified configuration, and compared with that of the previous configuration. If the error is minimized between the two configurations, then the orientations are retained, else the orientations are flipped to revert back to the previous configuration.