Multiscale Tomographic Wave–Matter Interaction Modeling to Enable Artifact-Free Material Defect Reconstruction

With only a few extraordinarily well-controlled exceptions, all materials contain defects [1]. Here, the term “defect” is used to describe localized regions featuring properties which vary from those of the surrounding material. Defects are commonly classified according to their morphology as voids, inclusions, cracks, etc. Such defects exist as ensembles spanning many length scales, from atomic-scale dislocations to macroscopic cracks and voids. An example of a defect created by lack of fusion in an additively manufactured part is illustrated in Fig. 1.

The presence of defects in engineered materials, components, and structures can have enormous deleterious performance ramifications [2]. As a result, an extensive variety of techniques have been developed to detect and identify defects in situ. Common approaches include ultrasonic [3–5], thermographic [6–8], shearographic [9–11], Lamb-wave [12–14], and X-ray tomographic (XCT) [15–17] methods. Similarly, a large body of work exists regarding the detection of defects in particular material classes, with some relevant areas of focus including metals [18–22], composites [23–28], electronics [29–32], and biologically derived materials [16,33–36].

The available methods for defect detection are based on waves interacting with materials bearing defects. These waves may be of a single wavelength (“monochromatic) or spectrally distributed (“polychromatic”). These waves interact via various mechanisms (refraction, scattering, etc.) with defects in a highly complex fashion. These interactions, combined with limitations in detector technology, create spurious, non-physical patterns or “artifacts” in the resulting scan data. The ultimate purpose of the present effort is to identify, implement, and validate a methodology for the removal of such artifacts from scan data, in order to accurately model and assess the performance impact of defects discovered in situ.

The present research is confined to the topic of XCT for the sake of practicality and conciseness. This choice was made due to the widespread application of XCT to performance-critical components and structures, particularly in the aerospace industry. However, the methods developed in the course of this research are based on a general form of the wave equation, and are readily extensible to other defect detection technologies. The remainder of this section is devoted to a description of XCT. This is followed by a motivational example, where the generation of scan artifacts in XCT data is demonstrated. It is shown that these artifacts introduce an unacceptable level of error in both the recovered defect geometry, and the resulting lifespan prediction based on that geometry and a simple cyclic fatigue model. In response to this shortcoming, a new method for modeling the interactions of defects with X-rays is presented. This method is based on the use of the Finite Element Method (FEM) to calculate solutions to the Helmholtz form of the wave equation. Combined with a novel defect parameterization method, it is shown that this modeling capability enables the calculation of artifact-free defect geometry utilizing a classic inversion method. A two-dimensional synthetic experimental dataset is used to demonstrate this capability. A road-map of the many future developments necessary to mature this preliminary research in to a broadly applicable methodology concludes the paper.

A significant part of the work in the present paper was originally presented at the 2020 ASME International Design Engineering Technology/Computers and Information in Engineering conferences [37].

The field of tomography originates with the mathematical work of Radon [38], which was first practically applied to the area of medical imaging in the late 1960s and early 1970s. Since that time XCT has become an invaluable tool in the medical diagnostic field, and has been adopted for a wide variety of additional applications. An excellent overview of XCT as applied to the areas of engineering and materials science is given by Maire and Withers [39]. A diagram of a generic XCT process is shown in Fig. 2. A wave emitter, in this case an X-ray source, generates an interrogating wave which is directed at the test object. Rotational and/or translational actuators position the specimen relative to the interrogation waves. After passing through the test object, the transmitted wave is measured by a (2D) detector, such as a scintillator crystal or a flat-panel sensor. This process is repeated for multiple orientations of the test object, most commonly a series of small rotations about a single axis. The individual transmission X-ray scans, combined with the pose information of the test object, are used by a reconstruction algorithm to produce an approximation of the three-dimensional object geometry.

From a computational perspective, the reconstruction algorithm is at the heart of XCT. The foundational method for XCT is the inverse radon transform, which is discussed in depth in Refs. [38,40]. While it remains the core of many XCT software tools, numerous refinements and alternatives have also been employed. Other common reconstruction algorithms include singular value decomposition [41], wavelet transforms [42], least-squares data fitting [43], and more recently, convolutional neural networks [44]. An interested reader may consult reference works such as Refs. [45,46] for a more detailed summary of available approaches.

XCT scans often exhibit noise and data artifacts [45]. For example, streaking in XCT scans is commonly observed, and is caused by combinations of undersampling (too few scans), photon starvation (inadequately energetic X-ray source), beam hardening (differential attenuation between different wavelengths of X-ray), and Compton scattering. Resultantly, many approaches have been developed to mitigate such deleterious phenomena. One common approach is filtered backprojection in the Radon transform, e.g., Refs. [47,48]. Additional approaches include those in Refs. [49,50]. The work of Jin et al. [49] is particularly relevant, as it demonstrates that incorporating even basic forward models of X-ray propagation permits scan reconstruction with greatly reduced artifacts.

There are several fundamental limitations associated with existing approaches for XCT filtering and noise-reduction limiting their applicability to the present work. The first of these limitations is due to the primarily bio-medical usage of XCT. The algorithms developed for filtering in this context typically smooth the output of the XCT data, in order to reduce visual noise and enable easier qualitative identification of scan features by human operators. Unfortunately, this process results in removing defect geometry information from the scan. Another issue common to existing techniques is the binarization of defect geometry into a pixel or voxel format. This step also removes important information from the scan data, particularly regarding small features. The section “Motivation” demonstrates these limitations on a simple test problem.

A notional example problem, demonstrating the need for improved, artifact-free defect reconstruction, is outlined in Fig. 3. A homogeneous circular domain (in 2D) is featured containing a single star-shaped defect. Synthetic tomographic data are generated by tracing rays from source points to corresponding detector points on opposing sides of the medium. The vector of detected intensities is denoted d. The circular domain containing the defect is rotated about its center, and successive scans are assembled into a matrix (S) which may be visualized as a sinogram as shown in Fig. 4. This sinogram clearly shows the trace of the defect’s presence.

The artifacts introduced by the reconstruction are immediately clear. The defect geometries are extracted via binarization at a 50% gray level threshold, and are shown in Fig. 6. These geometries were inserted in to a basic 2D finite-element model of a tensile test coupon in order to assess the impact of defect geometry errors on local stress/strain field intensification. The defect half-width is 100 μm, and defines a hole in a rectangular domain measuring 2000 μm wide by 3000 μm high. The domain size was chosen to be sufficiently large to approximate an infinite elastic plane. A unit-magnitude, uniform tensile load was applied at the top edge of the domain. A plane-stress formulation was used, and the domain was discretized into approximately 500,000 second-order triangular elements. A maximum mesh size of 0.1 μm was enforced in the near-field (a circle 150 μm radius centered on the defect). The material properties of mild structural steel were used, namely, E = 200 GPa and ν = 0.3. All conditions, except defect geometry, were kept constant between the four simulations. The mumps solver [52] was used to solve the resulting system of equations, using default parameters. The resulting distributions of strain energy density for each geometry are shown in Fig. 7.

The finite element model was also utilized to compute the fatigue life of a domain containing each of the four defect reconstructions. A S-N criterion was applied to the von-Mises stress with a Goodman mean stress correction. The differences in predicted fatigue life are given in Table 1.

From the comparisons of Figs. 6 and 7, and the data of Table 1, it is clear that these conventional approaches to tomographic reconstruction have the potential to introduce unacceptable errors in to the lifespan prediction. The best possible case overpredicts remaining life by nearly 50%, and more realistic scenarios produce predictions with an order of magnitude error. Other research, such as Ref. [53], has found similar levels of geometric error in XCT reconstructions used for industrial metrology. The use of contouring algorithms such as Ref. [54] may reduce the binarization artifacts appreciably, but other sources of error still remain.

The poor performance of the reconstruction in this respect is unsurprising for a number of reasons. First, it is noted that the overwhelming majority of XCT research has been carried out in the bio-medical context. In this context, the primary usage of XCT is diagnostic or qualitative in nature, and hence many of the derived approaches to tomograph reconstruction tend to over-smooth the underlying data in order to produce results that are easier for human interpretation. Secondarily, it is seen that the underlying pixel (or in 3D, voxel) representation of defect geometry is poorly suited to some specific defect shapes. For example, the extremely sharp cusps visible in Fig. 1 require an impractically high pixel resolution to capture stress intensification properly. Finally, it is clear that the heuristic filtering and smoothing operations commonly employed may only remove information, as they do not encapsulate any of the relevant physics of wave–defect interaction. From these observations it follows that a physics-aware tomographic reconstruction approach, which yields parametric representation of defect geometry, may offer a great benefit in the effort to quantitatively study the effects of defects on material/structural performance.

Motivated by the example given in the section “ Motivation”, a framework for tomographic reconstruction incorporating the physics of interactions between interrogating waves and defects is proposed. This methodology consists of three major components,

1. Mathematical parameterization of defects.

2. Forward modeling of the physical interactions of defects with interrogation waves.

3. Inversion of the forward model to extract defect geometry (in terms of the chosen parameterization) from experimental measurements.

The following sections describe each of these components in detail.

Prior work such as Ref. [49] utilizes the principles of ray optics, and ray tracing algorithms [58] to model the interaction of X-rays with materials in XCT. While such approaches are computationally expedient, consideration of the multi-scale geometric features of defects reveals a problem: many defects have features that are of a similar scale to the wavelength of interrogating X-rays. Examples of such features are the sharp cusps seen in Fig. 1. This necessitates the modeling of X-rays as wave phenomena, as illustrated conceptually in Fig. 9. Instead of rays propagating through various domains (Ω), the X-rays are modeled using the wave equation.

Equation (4) is discretized and solved using the FEM. Notably, this discretization may be performed in two or three dimensions, without reformulation of the underlying equations. For simplicity and computational expediency, the present work utilizes a 2D example problem and associated discretization. The development of more general and complete FE models in 3D, incorporating additional physical phenomena such as realistic source and detector behavior, will be a major area of effort in the future.

The section “Defect Interaction Modeling” defines a “forward” model. For a given combination of material, defect, and X-ray properties, the model predicts the X-ray intensity measured by the detector. For example, Fig. 11 shows a comparison between the detected intensity for an empty X-ray chamber, the same chamber containing a defect-free material, and a defect-bearing material. The large difference between the signals from non-defective and defective materials indicates that there is significant information about the defect structure conveyed in these data.

In order to quantitatively characterize the geometry of unknown defects from experimental X-ray data, this model must be inverted. This inversion may be accomplished via any of a large number of techniques [59]. In this case, a method previously developed by the authors for composite material characterization [60] is employed. First, the forward model is run i times, for multiple defect geometry parameter values, using a latin-hypercube sampling strategy. For each run, the detector response δ_i is recorded. Each such dataset consists of n_j discrete observations of field intensity (δ_ij). For the 2D example problem of the present work, these sampling points are taken along a line, corresponding to a 1D detector format. For the 3D case, the sampling points would be defined on a 2D plane, corresponding to the dimensions of a flat-panel or scintillating crystal X-ray detector.

The methodology outlined in the section “ Inversion Method” is generic in nature, and will require a great deal of refinement to become a practical and useful tool for defect analysis. In order to justify this effort, in this section the application of the methodology to a simple synthetic test problem is demonstrated.

The geometry illustrated in Fig. 10 was used, with a 30 mm circular material object bearing an off-center star-shaped defect. These domains were meshed using approximately 145,000 second-order triangular elements. As before, second-order triangular elements are employed. Elements within Ω_d are constrained to have an edge length no greater than 1 μm, those within Ω_m (which is a circle of 500 μm radius) are no larger than 10 μm, and those in Ω_e are no larger than 1000 mm. The Neumann boundary conditions of Eqs. (6) and (7) were enforced, with c_e = 343.4 and ρ_e = 1.2 (the values for air at standard temperature and pressure). For simplicity, the other parameters were defined c_m = 10c_e, c_d = 0.5c_e, ρ_m = 10ρ_e, and ρ_d = 0.5ρ_e. Any of the commonly available FEM solvers may be employed to solve the Helmholtz equation under these conditions. In the present work, the mathematica [61] computing system is employed for discretization and solution. Although the numerical solver used by mathematica is proprietary, it gives identical results to those calculated using the aforementioned sc mumps sparse/direct solver.

Solution of the FEM required 5 s on a laptop computer with a 6-core Intel Xeon processor (2.7 GHz clock frequency) and 64 Gb of system memory. Thee results of one run of the FEM are shown in Fig. 13. The defect parameters space were defined by 0.5 ≤ A ≤ 1.5, 0.5 ≤ B ≤ 1.5, 0.1 ≤ r ≤ 0.9, and r = s. This three-parameter space was sampled in a uniform, 5 × 5 × 5 grid to generate samples for the surrogate model. The surrogate model was defined by the linear interpolation of these data. A synthetic experimental dataset was generated by running the FEM model a final time, at A = B = 1.01, r = s = 0.013125 (a point not in the data used to generate the surrogate).

To solve the minimization problem the Nelder–Mead optimization algorithm [62] was employed. The optimizer’s reflection (α), expansion (γ), contraction (β), and shrinkage (ρ) coefficients were assigned the typical default values [63] of α = 1, γ = 2, β = 0.5, and ρ = 0.5. The optimization used 100 randomized starting points. The final result returned was by the optimization algorithm A = 1.01, B = 1.01, r = s = 0.0131336. This corresponds to an error of 0.06%. This result is two orders of magnitude better than that achieved using conventional approaches, as described previously. This very low error is a strong indication that the problem structure is well suited to inversion in this fashion, and suggests that further refinements of this approach are justified.

The previous sections have outlined a new method for the geometric reconstruction of defects from XCT data without the introduction of spurious artifacts. A motivational example was presented to illustrate the enormous effect that these artifacts have, resulting in especially severe over-estimation of structural fatigue life. To address this issue physics-based models, namely, solutions to the Helmholtz wave equation achieved using Finite Element methods, have been used to model the interactions of X-rays with materials bearing defects. This model utilizes a superquadratic parametric model of defect geometry. It has been shown that these techniques may be used to recover defect geometry from synthetic experimental data using a model inversion technique. The defect geometry recovered in this fashion exhibits no artifacts, and far lower error than is possible using existing methods based on ray tracing.

Although the present work has focused on XCT technology, the novel use of the general wave equation readily enables its use for other wave-based defect detection methods. Extension to such methods (i.e., ultrasonic mechanical waves) requires only the determination of the correct wave parameters (i.e., c and ρ) and the implementation of appropriate boundary conditions. The formulation employed by the present work is also readily applicable to three-dimensional defect reconstruction, requiring discretization and solution (but not reformulation) of the relevant wave equation in three dimensions. While the immediate work regarding the further development of this method will continue to focus on XCT, the generality required for extending to other techniques will be maintained and explored when possible.

The examples presented in this paper are extremely simple in nature, and necessarily omit many aspects of real XCT systems. In order to further develop the proposed method into a practical and useful tool, a great deal of future effort is required. Several of the most important modeling developments will include:

1. Enrichment of the wave model used to represent X-ray/material interactions. The modeling of various scattering phenomena will be particularly important.

2. Inclusion of more realistic source and detector physics including polychromatic wavelengths, and nonlinear transformations applied to the detected signal.

3. Extension of the model to three dimensions, with the goal of modeling cone-beam micro-tomography (uCT) systems commonly employed within the materials science context.

Further improvements to the other elements of the methodology will also be required. These may be conducted in parallel with the efforts listed previously. Important topics will include:

1. Identification of the best available defect geometry parameterizations. Possible approaches may include the superquadratic forms of the present work, B-rep geometry derived using NURBs or other geometry kernels, or volumetric fractal representations.

2. Selection of the best surrogate or reduced-order modeling approaches to employ for inversion. The linear interpolation of the present work will likely require excessively dense sampling when applied to geometries with larger numbers of parameters.

3. Exploration of sequential or adaptive sampling schemes for generating the data for building the surrogate models.

4. Selection of the best numerical optimization tools for solving the minimization problem necessary to recover the defect geometry.

In addition to these extensions and improvements, rigorous experimental validation of the methodology will be conducted. This will be accomplished by the fabrication of test specimens containing defects of known geometry. Composite coupons, likely of the well-characterized AS4/3501-6 carbon-epoxy system, containing inclusions of polymer which are highly translucent to X-ray radiation will be manufactured. The geometry of these inclusions will be precisely measured using appropriate 3D-scanning or microscopic metrology techniques before specimen fabrication. The specimens will then be scanned using a highly accurate Zeiss Xradia 520 Versa micro-tomography system (uCT). The data gathered from these scans will be analyzed via the inversion methodology outlined in this work. Validation will be achieved by comparison of the predicted geometry to the actual defect geometry. Importantly, these predictions will be made in a “blind” fashion in order to prevent over-fitting of the surrogate models. Multiple validation experiments containing representative defect geometries will demonstrate that the methodology is sufficiently general. Additionally, after the blind comparison has been made the test articles will be measured using a destructive technique such as serial sectioning. This will provide a second independent validation that the methodology predicts defect geometry with sufficiently high accuracy. Additionally, this step will ensure that no unintentional defects have been introduced into the validation experiment.

Despite the magnitude of these future efforts, the proposed methodology offers unique and potentially important features. The ability to extract artifact-free, high accuracy defect geometries from XCT and uCT data will offer great utility in the lifespan estimation of components and structures. Importantly, this will not require extensive “human in the loop” scan segmentation efforts, as is required by ongoing machine-learning research. Therefore this capability may also be utilized in high-throughput experimental contexts, such as automated acceptance and qualification testing. With additional effort, the methodology may eventually become broadly applicable to many wave-based non-destructive testing technologies.

The authors acknowledge support for this work by the Office of Naval Research through the Naval Research Laboratory’s core funding.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.