ELAS3D-Xtal: An OpenMP-accelerated crystal elasticity solver with automated experiment-driven microstructure generation

TL;DR

ELAS3D-Xtal delivers a high-performance, OpenMP-accelerated solver for 3D elasticity in defect-containing polycrystals, combining a matrix-free PCG method with a block-Jacobi preconditioner and rotated stiffness tensors to handle crystal anisotropy. The framework includes an integrated native microstructure generator, voxel-to-grain assignment, and pore incorporation from XCT data, with HDF5 outputs for scalable post-processing. Validation against the Eshelby inclusion problem confirms accuracy, while application to LPBF SS316L demonstrates how elastic anisotropy and defect morphology govern local stress concentrations. The work enables rapid, high-resolution defect-mechanics studies on commodity hardware and sets the stage for future multi-fidelity fatigue assessments linking processing-induced defects to reliability.

Abstract

This paper introduces ELAS3D-Xtal, a high-performance Fortran/OpenMP upgrade of the NIST ELAS3D voxel-based finite element solver for computing 3D elastic fields in polycrystals with defects. The code supports crystal anisotropy by precomputing rotated stiffness tensors from user-specified orientations and solves the equilibrium problem with a matrix-free, OpenMP-parallel preconditioned conjugate-gradient (PCG) method using a point-block Jacobi preconditioner. On a single shared-memory multicore PC, OpenMP threading accelerates the baseline CG solver by ~10X, while the block-preconditioned CG solver achieves 53-61X speedup relative to the serial CG baseline for meshes from 100^3 to 500^3 voxels (scaling to domains up to 800^3 voxels). Accuracy is validated against the analytical Eshelby inclusion solution. ELAS3D-Xtal also integrates microstructure construction, including statistically calibrated polycrystal generation via spatial filtering and parallel voxel-to-grain assignment, direct pore insertion from XCT centroid/radius data, and texture assignment. Full-field phase, orientation, and stress outputs are written in HDF5 to enable scalable post-processing and defect-mechanics workflows. Applications are demonstrated for (i) anisotropy-controlled defect-scale stress fields and (ii) LPBF SS316L microstructures with gas, lack-of-fusion, and keyhole pore morphologies.

Figures (16)

• Figure 1: Computational workflow for ELAS3D-Xtal. The algorithm proceeds as follows: (1--2) Initialization & Prep: Simulation parameters are defined and microstructure data is read from HDF5 inputs to build periodic neighbor tables. (3) Stiffness Assembly: A decision tree directs OpenMP-parallelized assembly for either isotropic media (flag=0, uniform moduli) or polycrystals (flag=1, rotated stiffness tensors). (4--5) Solver Loop: After applying macroscopic strain and initializing energy functionals, the OpenMP-accelerated PCG solver employs point-block Jacobi preconditioning and Fletcher-Reeves updates until convergence. (6) Output: Final phase maps and full-field stresses are exported.
• Figure 2: Computational workflow for microstructure generation. The algorithm integrates experimental pore data from XCT scans with statistically representative grain morphologies. Spatial filtering and parallel block processing enable efficient voxel-to-grain assignment. Pore integration overwrites grain assignments with spherical defect regions. Crystallographic orientations are assigned using either fiber texture or random sampling.
• Figure 3: Schematic of the microstructure reconstruction framework. The pipeline features: (1) direct integration of pore data from XCT scans; (2) efficient parallel processing via spatial filtering; and (3) flexible orientation assignment allowing for both uniform SO(3) sampling and Gaussian-spread fiber textures.
• Figure 4: Simulation setup for the 3D Eshelby inclusion problem under macroscopic strain loading along the $X_3$-direction. Material properties are defined for an elastic matrix ($E_m=194$ GPa, $\nu_m = 0.2934$, colored white) surrounding a zero-stiffness void inclusion ($E_{in}=0.0$ GPa, $\nu_{in} = 0.0$, colored blue). (a) 3D phase map visualization of the domain. (b) Voxelized representation utilizing filled cubes to emphasize the discrete computational grid. (c) 2D cross-sectional view taken at the central $X_1$-plane, illustrating the phase distribution overlaid with the finite element mesh grid.
• Figure 5: Validation of the numerical solver against the analytical Eshelby solution using a high-resolution $500^3$ mesh. The columns display the analytical solution (left), the numerical PCG solution (middle), and the difference field (right). The rows correspond to central cross-sections in three orthogonal planes: (a) the $X_2-X_3$ plane at the center of $X_1$; (b) the $X_1-X_2$ plane at the center of $X_3$; and (c) the $X_1-X_3$ plane at the center of $X_2$.