Neutron Diffraction Strain Tomography: Demonstration and Proof-of-Concept

TL;DR

This paper explores the feasibility of neutron diffraction strain tomography from diffraction measurements, developing a reconstruction algorithm and demonstrating its effectiveness through simulation and experimental validation.

Contribution

It introduces a novel method for strain tomography using diffraction data and provides a proof-of-concept demonstration with validation against traditional methods.

Findings

Successful reconstruction of strain fields in simulations

Experimental validation on an axisymmetric sample

Algorithm shows promise for non-destructive strain analysis

Abstract

Recently, a number of reconstruction algorithms have been presented for residual strain tomography from Bragg-edge neutron transmission measurements. In this paper, we examine whether strain tomography can also be achieved from diffraction measurements. We outline the proposed method and develop a suitable reconstruction algorithm. This technique is demonstrated in simulation and a proof-of-concept experiment is carried out where the strain field in an axisymmetric sample is reconstructed and validated against conventional diffraction strain scans.

Equations

$$\lambda =2d\sin \vartheta ,$$

$$⟨ϵ⟩=\frac{d-d_0}{d_0},$$

$$\begin{gathered} ⟨ϵ⟩=\frac{1}{{\int }_0^L{e}^{-\mu L_T(s)}\text{d}s}{\int }_0^L{e}^{-\mu L_T(s)}{\hat{\mathbf{\kappa }}}^{\mathsf{T}}\mathbf{ϵ}(s) \\ \hat{\mathbf{\kappa }}\text{d}s. \end{gathered}$$

$$\mathbf{m}(\mathbf{x})$$

$$\mathbf{K}(\mathbf{x},{\mathbf{x}}^{\prime })$$

$$\mathbf{K}(\mathbf{x},{\mathbf{x}}^{\prime })={\sigma }_f^2\text{exp}\left(\frac{-{\Vert \mathbf{x}-{\mathbf{x}}^{\prime }\Vert }^2}{2{\ell }^2}\right).$$

$$y_i={\mathcal{L}}_{{\mathbf{\eta }}_i}\mathbf{f}(\mathbf{x})+e_i.$$

$$\left[\begin{array}{c}\mathbf{Y}\\ \hat{\mathbf{f}}({\mathbf{x}}_{\ast })\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mathbf{\mu }}_{\mathbf{y}}\\ \mathbf{m}({\mathbf{x}}_{\ast })\end{array}\right],\left[\begin{array}{cc}{\mathbf{K}}_{\mathbf{y}{\mathbf{y}}^{\prime }}+{\mathbf{\Sigma }}_{\mathbf{m}}& {\mathbf{K}}_{\mathbf{y}{\hat{\mathbf{f}}}^{\prime }}\\ {\mathbf{K}}_{\hat{\mathbf{f}}{\mathbf{y}}^{\prime }}& \mathbf{K}({\mathbf{x}}_{\ast },{\mathbf{x}}_{\ast })\end{array}\right]\right)$$

$${\mathbf{K}}_{\mathbf{y}{\hat{\mathbf{f}}}^{\prime }}=\left[\begin{array}{c}{\mathcal{L}}_{{\mathbf{\eta }}_1}\mathbf{K}(\mathbf{x},{\mathbf{x}}_{\ast })\\ ⋮\\ {\mathcal{L}}_{{\mathbf{\eta }}_n}\mathbf{K}(\mathbf{x},{\mathbf{x}}_{\ast })\end{array}\right]$$

$${\mathbf{K}}_{\mathbf{y}{\mathbf{y}}^{\prime }}=\left[\begin{array}{ccc}{\mathcal{L}}_{{\mathbf{\eta }}_1}\mathbf{K}(\mathbf{x},{\mathbf{x}}_{\ast }){\mathcal{L}}_{{\mathbf{\eta }}_1}^{\prime \mathsf{T}}& \cdots & {\mathcal{L}}_{{\mathbf{\eta }}_1}\mathbf{K}(\mathbf{x},{\mathbf{x}}_{\ast }){\mathcal{L}}_{{\mathbf{\eta }}_n}^{\prime \mathsf{T}}\\ ⋮& \ddots & ⋮\\ {\mathcal{L}}_{{\mathbf{\eta }}_n}\mathbf{K}(\mathbf{x},{\mathbf{x}}_{\ast }){\mathcal{L}}_{{\mathbf{\eta }}_1}^{\prime \mathsf{T}}& \cdots & {\mathcal{L}}_{{\mathbf{\eta }}_n}\mathbf{K}(\mathbf{x},{\mathbf{x}}_{\ast }){\mathcal{L}}_{{\mathbf{\eta }}_n}^{\prime \mathsf{T}}\end{array}\right].$$

$${\mathbf{\mu }}_{{\mathbf{f}}_{\ast }|\mathbf{Y}}$$

$${\Sigma }_{{\mathbf{f}}_{\ast }|\mathbf{Y}}$$

$$\frac{\text{d}\sigma_{rr}}{\text{d}r}+\frac{1}{r}\bigg{(}\sigma_{rr}(r)-\sigma_{\theta\theta}(r)\bigg{)}=0.$$

$$\mathbf{\sigma }(r)=\left[\begin{array}{c}{\sigma }_{rr}(r)\\ {\sigma }_{\theta \theta }(r)\end{array}\right]=\left[\begin{array}{c}\varphi (r)\\ r\frac{\text{d}\varphi }{\text{d}r}+\varphi (r)\end{array}\right]$$

$$\mathbf{ϵ}(r)=\left[\begin{array}{c}{ϵ}_{rr}(r)\\ {ϵ}_{\theta \theta }(r)\end{array}\right]=\frac{1}{E}\left[\begin{array}{c}{\sigma }_{rr}(r)-\nu {\sigma }_{\theta \theta }(r)\\ -\nu {\sigma }_{rr}(r)+{\sigma }_{\theta \theta }(r)\end{array}\right]$$

$$\mathbf{ϵ}(r)=\frac{1}{E}\left[\begin{array}{c}(1-\nu )-\nu r\frac{\text{d}}{\text{d}r}\\ (1-\nu )+r\frac{\text{d}}{\text{d}r}\end{array}\right]\varphi (r).$$

$$\mathbf{ϵ}(r)=\left[\begin{array}{c}(1-\nu )-\nu r\frac{\text{d}}{\text{d}r}\\ (1-\nu )+r\frac{\text{d}}{\text{d}r}\end{array}\right]\varphi (r)=\mathbf{H}\varphi (r).$$

$$⟨ϵ⟩={\left[\begin{array}{c}{W}_{\mu }{\int }_0^L{w}_{\mu }(s)w_r(s)(\cdot )\text{d}s\\ W_{\mu }{\int }_0^L{w}_{\mu }(s)w_{\theta }(s)(\cdot )\text{d}s\end{array}\right]}^{\mathsf{T}}\mathbf{ϵ}(r)=\mathbf{M}\mathbf{ϵ}(r).$$

$$⟨ϵ⟩=\mathbf{M}\mathbf{ϵ}=\mathbf{MH}\varphi (r)={\mathcal{L}}_{\eta }\varphi (r)$$

$$\underset{\boldsymbol{h}}{\operatorname{arg}\,\operatorname{max}}\;\Big{(}\log\text{det}(\boldsymbol{K}_{\boldsymbol{yy}}(\boldsymbol{h})+\boldsymbol{\Sigma}_{m})-\boldsymbol{y}^{\mathsf{T}}(\boldsymbol{K}_{\boldsymbol{yy}}(\boldsymbol{h})+\boldsymbol{\Sigma}_{m})^{-1}\boldsymbol{y}\Big{)}.$$

$$\left[\begin{array}{c}{ϵ}_{rr}(r)\\ {ϵ}_{\theta \theta }(r)\end{array}\right]=e_0\left[\begin{array}{c}\frac{(7+5\nu )R^2+(1+\nu )(9r^2-16Rr)}{12R^2}-{\left(1-\frac{r}{R}\right)}^2\\ \frac{(7+5\nu )R^2+(1+\nu )\left(3r^2-8Rr\right)}{12R^2}-{\left(1-\frac{r}{R}\right)}^2\end{array}\right],$$

$$⟨ϵ⟩=\frac{1}{{\int }_0^L{e}^{-\mu L_T(s)}\text{d}s}{\int }_0^L{e}^{-\mu L_T(s)}{\hat{\mathbf{\kappa }}}^{\mathsf{T}}\mathbf{ϵ}(s)\hat{\mathbf{\kappa }}\text{d}s,$$

$${\hat{\mathbf{\kappa }}}^{\mathsf{T}}\mathbf{ϵ}(s)\hat{\mathbf{\kappa }}={\kappa }_x^2{ϵ}_{xx}+2{\kappa }_x{\kappa }_y{ϵ}_{xy}+{\kappa }_y^2{ϵ}_{yy}.$$

$$\hat{\mathbf{\kappa }}=\left[\begin{array}{cc}\cos (\pi /2-\vartheta )& \sin (\pi /2-\vartheta )\\ -\sin (\pi /2-\vartheta )& \cos (\pi /2-\vartheta )\end{array}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]=\left[\begin{array}{c}\cos \vartheta \\ \sin \vartheta \end{array}\right].$$

$${\hat{\mathbf{\kappa }}}^{\mathsf{T}}\mathbf{ϵ}(s)\hat{\mathbf{\kappa }}$$

$$\displaystyle\qquad+2\sin\vartheta\cos\vartheta\bigg{(}\sin{\theta}\cos{\theta}(\epsilon_{rr}(r)-\epsilon_{\theta\theta}(r))\bigg{)}$$

$$\displaystyle\qquad\qquad+\sin^{2}\vartheta\bigg{(}\epsilon_{rr}(r)\sin^{2}{\theta}+\epsilon_{\theta\theta}(r)\cos^{2}{\theta}\bigg{)}$$

$$={ϵ}_{rr}(r){\cos }^2(\vartheta -\theta )+{ϵ}_{\theta \theta }(r){\sin }^2(\vartheta -\theta )$$

$$\langle\epsilon\rangle\approx W_{\mu}\int_{0}^{L}w_{\mu}(s)\bigg{(}w_{r}(s)\epsilon_{rr}(r)+w_{\theta}(s)\epsilon_{\theta\theta}(r)\bigg{)}\;\text{d}s$$

$$r=r(s)=\sqrt{x^2+(s-\sqrt{2pR-p^2}{)}^2}$$

$$w_{\mu }(s)$$

$$W_{\mu }$$

$$w_r(s)$$

$$w_{\theta }(s)$$

Full text

††journal:

1 Introduction

For many years, neutron diffraction instruments have been able to provide high-precision measurements of strain (down to $\sigma=1\times 10^{-5}$ uncertainty). by observing changes in the atomic lattice spacing $d$ within a polycrystalline samples [1]. These measurements can be performed with resolution down to $0.5$mm3 based on practical limitations in gauge volume size.

These measurements rely on Bragg’s law, which provides the condition for constructive interference of neutron radiation in a lattice:

[TABLE]

where $\lambda$ is the neutron wavelength, $d$ the average spacing of all planes within the gauge volume aligned with the direction $\boldsymbol{\hat{\kappa}}$ bisecting the incident and diffracted beams, and $\vartheta$ the half-angle of diffraction, as shown in Figure 1.

In conventional neutron strain scanning, $\lambda$ is known and $\vartheta$ measured or vice-versa and the average strain in the $\boldsymbol{\hat{\kappa}}$ direction is then calculated by:

[TABLE]

where $d_{0}$ is the lattice spacing in a reference (unstrained) sample.

By taking a number of measurements from different sample orientations, it is possible to resolve the entire triaxial strain tensor at one or a number of points within a body. This ‘pointwise’ approach has obvious application where the area of interest is known; e.g. probing the stress at a fixed location as a function of some external stimulus, or examining the effects of a crack/weld/other feature.

Where full-field strain maps are required (e.g. when evaluating the residual strains locked-in by novel manufacturing techniques), it is commonplace to interpolate between sets of discrete, distributed, point-wise measurements. Such interpolations present two significant issues:

The quality of the interpolation is dependent on the resolution of the pointwise measurements. Achieving sufficient resolution (particularly in 3 dimensions) is not always practical as beamtime is a limited resource. 2. 2.

Interpolations are not guaranteed to represent the underlying field with a high level of fidelity — i.e. they may not satisfy physical constraints such as equilibrium or compatibility.

In recent years, a number of algorithms for full-field strain tomography from neutron transmission measurements have been developed and have shown promise both in simulation and on experimental data [3, 4, 5, 6, 7, 8, 9, 10, 11, 12].

These algorithms rely on Bragg-edge imaging, whereby the the transmitted (as opposed to diffracted) neutrons are counted [13]. Presently, only a handful of facilities are suitable for performing strain tomography using these techniques.

We propose that strain imaging is also possible using diffraction geometry and instrumentation, and that the recently developed reconstruction algorithms can be extended to this problem.

If successful, this will significantly extend the utility of these algorithms in two major ways. First, this may allow strain tomography at the multitude of existing diffraction strain scanners around the world, and secondly, this may open the door for a combined reconstruction algorithm which utilises transmission and diffraction measurements in tandem. This tantalising prospect has the potential to reduce beamtime requirements and improve the quality of strain reconstructions.

In this paper, we present a method and algorithm for strain tomography using diffraction geometry. As a proof-of-concept, we restrict this demonstration to axisymmetric systems. We first demonstrate the approach using simulated measurements, and then reconstruct the strain within a crushed, axisymmetric disc from experimental data. We conclude by briefly discussing potential improvements to the algorithm and the considerations in extending this method to arbitrary 2D systems.

2 Method

2.1 Experimental Setup

We propose that conventional diffraction strain scanners can be utilised to obtain ‘ray-like’ measurements and perform strain imaging analogous to a transmission setup. The forward mapping has been examined in [14]. This technique may be feasible for both constant-wavelength (e.g. see Figure 2) and time-of-flight strain scanners. While the finer details will differ, we provide a general overview for the former:

A polychromatic beam of neutrons is generated by a reactor source and is directed toward the instrument via a shielded beam line. 2. 2.

This beam is simultaneously redirected, focussed, and reduced to a single wavelength using a curved monochromator. 3. 3.

The monochromatic beam floods the sample (i.e. a fully-open primary slit) and neutrons are scattered outward in various directions by the many, differently oriented planes according to Bragg’s law. 4. 4.

A cadmium secondary slit defines an effective gauge volume in the form of a long, thin region that can be approximated as a ray111Obviously, the secondary slit width $w$ has an effect on the validity of the ray assumption. This is discussed in Section 4. — see Figure 3. 5. 5.

Analogous to conventional strain scanning, the measured lattice spacing $d$ then provides the average normal strain in the $\boldsymbol{\hat{\kappa}}$ direction along this ray.

By sweeping the sample past the secondary slit (essentially varying $x$ in Figure 3), it is possible to obtain a profile of measurements analogous to those which can be measured using pixelated detectors in a transmission setup.222For time-of-flight instruments, it may be possible to achieve something equivalent using a collimator and imaging detector.

Many of the algorithms already developed for Bragg-edge transmission tomography may be adapted to this geometry. Two significant differences exist;

In this regime, the measurement direction is not aligned with the direction of a ray. 2. 2.

The incident and diffracted beams are not collinear, meaning attenuation along the neutron ray must be considered.

Attenuation provides a significant challenge — each measurement now represents a weighted average of the strains along a ray. The contribution from each point along the ray is dependent on the total path length seen by a neutron diffracted at that point.

2.2 Measurement Model

Considering the effect of attenuation, we propose that diffraction tomography measurements can be modelled by a ray transform of the form:

[TABLE]

For a given measurement, strain in the $\boldsymbol{\hat{\kappa}}$ direction is averaged along the ray with length $L$ as shown in Figure 3. $L_{T}(s)$ is the total path length seen by a neutron diffracted at a point $s$ units along this ray (including the incident path). $\mu$ is the effective neutron attenuation coefficient of the material being measured.

2.3 Algorithm Selection

A number of algorithms [3, 4, 5, 6, 7, 8, 9, 10, 11, 12] have recently been presented for transmission strain tomography, and most could be modified to suit the problem at hand.

Of the those presented thus far, the most promising algorithms have been those which model the unknown strain field as a Gaussian Process (GP). Compared to other approaches, the convergence of GP regression algorithms for strain reconstruction has generally been significantly more rapid, partially owing to their implicit implementation of the equilibrium constraint and explicit implementation of known boundary (loading) conditions [11].

GP regression also presents a number of other benefits compared to other algorithms. Being a probabilistic approach, GPs are able to utilise the confidence of each measurement in the resulting reconstruction (i.e. by giving less weighting to more uncertain measurements). The GP technique also allows the confidence interval of the reconstruction to be calculated at every point.

Finally, the GP technique is non-parametric. This means that complexity does not increase with the resolution of the reconstruction. While a small, fixed number of so-called hyperparameters influence the reconstruction, these can be optimised from the measurements alone using a likelihood maximisation method.

For these reasons, we will implement a GP reconstruction algorithm in this paper.

2.4 Brief Introduction to Gaussian Processes

GP regression is explained in detail in [15]. Specifics related to implementing this technique for strain reconstruction can be found in [11, 12].

Briefly, a GP is a machine learning technique that models an unknown field as a Gaussian random distribution of functions $\boldsymbol{f(x)},\;\boldsymbol{x}\in\mathbb{R}^{\text{dim}(\boldsymbol{x})}$. This distribution is characterised by a mean function $\boldsymbol{m(x)}$ and covariance function $\boldsymbol{K(x,x^{\prime})}$, where;

[TABLE]

The choice of covariance function $\boldsymbol{K(x,x^{\prime})}$ can have a significant impact on the resulting reconstruction. For this application, we build a GP around the squared-exponential covariance function, which assumes a high degree of smoothness and has shown promise in modelling strain [16, 11, 12] — a usually smooth phenomena;

[TABLE]

Here, $\sigma_{f}^{2}$ is the variance on the prior (sometimes called the signal variance), and $\ell$ is a length-scale. These hyperparameters inform the most likely functions to be drawn from the GP. For instance, a small length-scale favours rapidly changing reconstructions, while larger length-scales motivate smoother, slower-varying realisations.

Gaussian process regression involves estimating a function value at a user specified sample point $\boldsymbol{x_{*}}$ given a set of data, $\mathcal{D}=\left\{y_{i},\boldsymbol{\eta}_{i}\;|\;\forall\;i=1,\dots,n\right\}$, where each measurement is of the form;

[TABLE]

Here, $\boldsymbol{\mathcal{L}}_{\boldsymbol{\eta}_{i}}\boldsymbol{f(x)}$ is a linear transformation of the underlying function $\boldsymbol{f(x)}$ that is parametrised by the set $\boldsymbol{\eta}_{i}$. $e_{i}\sim\mathcal{N}(0,\sigma_{i}^{2})$ is assumed zero-mean Gaussian noise with standard deviation $\sigma_{i}$.

Given that GPs are closed under linear operators [17, 18], the measurements $\boldsymbol{Y}=[y_{1},y_{2},\ldots y_{n}]^{\mathsf{T}}$ and a function value estimate $\hat{f}(\boldsymbol{x_{*}})$ are jointly Gaussian [15];

[TABLE]

Where $\boldsymbol{\Sigma_{m}}$ is a diagonal matrix in which the $i^{\text{th}}$ entry corresponds to the variance of the $i^{\text{th}}$ measurement, and the cross-covariance matrix $\boldsymbol{K_{y\hat{f}^{\prime}}}$ and covariance matrix $\boldsymbol{K_{yy^{\prime}}}$ are given by333note: $\boldsymbol{K_{\hat{f}y^{\prime}}}=\boldsymbol{K_{y\hat{f}^{\prime}}}^{\mathsf{T}}$:

[TABLE]

and

[TABLE]

Here, the notation $\boldsymbol{\mathcal{L}}^{\prime}$ is used to distinguish a linear transform operating on $\boldsymbol{f(x^{\prime})}$ from a linear transform $\boldsymbol{\mathcal{L}}$ operating on $\boldsymbol{f(x)}$. Finally, we can condition the prior estimate $f(\boldsymbol{x_{*}})$ on the known measurement values to give a posterior estimate with mean and variance obtained from;

[TABLE]

2.5 A Gaussian Process for Diffraction Tomography

In this paper, we develop an algorithm for axisymmetric, 2D systems. This choice is motivated by the simplifications afforded by geometry and the reduced measurement time required for validation in this proof-of-concept study.

Conceptually, there is well defined path to extend this method to arbitrary 2D and even 3D systems. Practical issues may arise — these are discussed in Section 4.

In the process, we also contribute a GP for strain which leverages the strong constraint provided by axisymmetry — this formulation may be useful outside this specific problem.

We assume a 2D circular sample of radius $R$ with origin at it’s centre and are concerned with reconstructing the tensor strain distribution, $\boldsymbol{\epsilon}$, within this sample. Naturally, we will work in polar coordinates $(r,\theta)$. Under an axisymmetry assumption, the strain tensor can be written in terms of two nonzero in-plane components which vary only in the radial direction: $\boldsymbol{\epsilon}(r)=\begin{bmatrix}\epsilon_{rr}(r)&\epsilon_{\theta\theta}(r)\end{bmatrix}^{\mathsf{T}}.$

Our problem is constrained by equilibrium, which is described in terms of the stress tensor, $\boldsymbol{\sigma}(r)=\begin{bmatrix}\sigma_{rr}(r)&\sigma_{\theta\theta}(r)\end{bmatrix}^{\mathsf{T}}$ — a linear transformation of strain — in polar coordinates by the differential equation:

[TABLE]

To encode this constraint, we define stress via a scalar potential $\phi(r)$ through the following mapping:

[TABLE]

By relating the components of stress in this way, any stress field resulting from the potential $\phi(r)$ automatically satisfies equilibrium444Direct substitution into Equation 2 confirms this..

This approach closely resembles the Airy stress function technique first encoded in a GP in [12], but avoids the singularity near $r=0$ present in the Airy stress mapping for polar coordinates. The mapping chosen here also encodes a key constraint that arises from the axisymmetric polar coordinate system; $\sigma_{rr}(0)=\sigma_{\theta\theta}(0)$.

In this paper, we develop our algorithm assuming plane-stress555A minor modification is required for plane-strain., for which Hooke’s Law takes the form:

[TABLE]

We then have the following relationships between the scalar potential for which we build our GP, $\phi(r)$, and the strains we wish to estimate:

[TABLE]

To improve numerical stability, we neglect the scaling factor $1/E$ and relate the potential to the strains through the mapping $\boldsymbol{H}$ by:

[TABLE]

The strain field is mapped to a measurement through the measurement model:

[TABLE]

With weights $W_{\mu},w_{\mu}(s),w_{r}(s)$ and $w_{\theta}(s)$ as defined in Appendix A. These mappings $\boldsymbol{H}$ and $\boldsymbol{M}$ together comprise the linear transformation $\boldsymbol{\mathcal{L}}_{\eta}$ that acts on $\phi(r)$ to provide a measurement $\langle\epsilon\rangle$:

[TABLE]

2.6 Hyperparameter Selection

The GP hyperparameters $\boldsymbol{h}=\begin{bmatrix}\ell&\sigma_{f}\end{bmatrix}^{\mathsf{T}}$ can have a significant impact on the fidelity of a reconstruction. Selection of these hyperparameters does not require a-priori knowledge of the system, or user intervention — they can be determined from the measurements alone.

This is achieved by posing an optimisation problem to maximise the marginal likelihood of the measurements given a set of hyperparameters $\boldsymbol{h}$ [15]. In practice, for numerical stability, the hyperparameters are estimated by maximising the log marginal likelihood:

[TABLE]

There may be several local minima for a given set of measurements. These can correspond to different interpretations of the data. For instance, the same data may involve noisy measurements of a smooth function or precise measurements of a quickly-varying function. Use of a multi-start optimisation procedure and/or non-gradient based method (e.g. simulated annealing) can help to avoid these local minima and find the hyperparameters that are most likely given a set of measurements.

3 Demonstration

3.1 Simulation

We first demonstrate our algorithm for the axisymmetric strain field previously examined in [9]:

[TABLE]

where $e_{0}$ is a scaling factor and $R$ the sample radius.

This field satisfies equilibrium but not compatibility (i.e. is residual), and corresponds to the elastic component of strain set up in response to a hydrostatic eigenstrain of the form $\epsilon^{*}_{rr}(r)=\epsilon^{*}_{\theta\theta}(r)=e_{0}(1-\tfrac{r}{R})^{2}.$ The method by which this field was generated is described in detain in the appendix of [9].

A total of 15 diffraction tomography measurements were simulated across the width of a sample with radius $R=6.5$mm via Equation 1 to form the profile shown in Figure 5(a). Material properties and attenutation characteristics of steel were chosen: $E=220$ GPa, $\nu=0.3$, $\mu=120$ m*-1*. Measurements were corrupted by simulated mean-zero Gaussian measurement noise with standard deviation $\sigma=0.5\times 10^{-4}$ — typical of that which can be expected on a diffraction strain scanner given gauge volumes of this size and sensible measurement times [19].

Hyperparameters $\boldsymbol{h}=\begin{bmatrix}\ell&\sigma_{f}\end{bmatrix}^{\mathsf{T}}$ of the covariance function were determined by the likelihood maximisation procedure described earlier and were found to be 11.2mm and $6.7\times 10^{-3}$, respectively.

Reconstruction results are shown in Figure 5(b) and (c), and show close agreement with the true field, which lies within the 1-$\sigma$ confidence interval of the reconstruction over it’s entire span.

The reconstruction differs most significantly from the true field near the sample centre (at small $r$). The confidence interval of the result is also widest here. The cause seems to be twofold:

Strains near $r=0$ do not contribute at all to the majority of measurements as only rays passing close to the sample centre can see these strains (the smallest $r$ seen by a given ray is $|x|$). 2. 2.

Those rays which do manage to probe strains at small $r$ do so with large path lengths and consequently the relative weighting of points near $r=0$ in the average is small.

The reconstruction was seen to converge to the true field as additional measurements were added near $x=0$ and as less attenuation was simulated.

Overall, this simulation shows a similar quality reconstruction as [9], however it should be noted that this algorithm required far fewer measurements (15 vs 512) to achieve this result. This is largely owed to the intrinsic equilibrium constraint encoded in the GP.

3.2 Experiment

Following success in simulation, the reconstruction algorithm was applied to data collected on the KOWARI constant wavelength diffractometer within the Australian Centre for Neutron Scattering (ACNS) at the Australian Nuclear Science and Technology Organisation (ANSTO) [19, 20].

Sample Design

The sample consisted of a stack of small aluminium discs, which, with respect to Figure 5, were manufactured as follows.

1. (a)

Five stepped discs were manufactured from $D\approx 12$mm diameter 7075-T6 aluminium round bar with dimensions as shown. 2. (b)

These discs were then individually crushed in a deformation-controlled process to a uniform thickness of $2.89\pm 0.02$ mm. 3. (c)

The crushed discs were then each symmetrically faced to a thickness of $t=1$ mm. 4. (d)

Finally, the discs were stacked and glued together to form one sample.

The final dimensions of each disc are in line with the typical rule-of-thumb for plane stress systems ($D/t\approx 10$) [21]. Stacking was essential to provide sufficient material for both validation and diffraction tomography measurements to be performed in the limited available beamtime.

Calibration and Experimental Setup

Initial work centred on characterising the beam and finding an optimal focussing condition for the monochromator on KOWARI.

Intensity variation across the beam introduces an additional non-uniform weighting to the measured average strain along a ray. The intensity profile was characterised using an iron-powder standard sample and tracking the Fe (110) reflection. The position of this peak is similar to that of the Al (200) reflection used for the tomography measurements, however Fe provided faster measurements.

By adjusting the monochromator focus, variation in this profile was minimised. The final condition achieved less than 10% variation across the sample and was assumed constant.

Nominally 90-degree geometry was adopted (89.6∘), and after line scans to find the sample edges, 11 diffraction tomography measurements were taken across the width of the sample using a $w=1$ mm secondary slit width ($1\times L\times 5$ mm gauge volumes). Measurements were performed until 5000 counts had been recorded at each point. Total measurement time was approximately 12 hours.

Validation measurements of the hoop, radial and axial components of strain were performed using a $1\times 1\times 5$ mm gauge volume and tracking the Al (311) reflection and were also performed over 12 hours. Results from the two experiments were correlated using a standard powder sample which was measured under both regimes. The unstrained lattice parameter $d_{0}$ was calculated from the three measured strain components using the plane-stress assumption.

Results

Reconstruction results are shown in Figure 6. In general, the reconstruction shows excellent agreement with the validation measurements — with almost all measurements lying within the 1-$\sigma$ confidence interval of the reconstruction.

Quantitatively, the difference between the tomographic reconstruction and validation measurements was Gaussian, with mean $\approx-5.5\times 10^{-5}$ and standard deviation $\approx 2.8\times 10^{-4}$ — respectively two and one order of magnitude less than the strains under examination. This suggests that the difference is largely due to unrejected measurement noise and not a systematic bias. Note that the systematic error near $r=0$ due to attenuation seen in the simulated example are not present here — aluminium has an attenuation coefficient nearly ten times smaller than steel.

4 Discussion and Potential Future Extensions

This proof-of-concept demonstration leaves many potential improvements to the technique on the table for future implementation:

A sensitivity analysis was performed as to the effect of secondary-slit width on simulated tomography profiles. It was found that the variation was minimal up to slit widths of 2 mm for some smooth strain fields. However, this effect may be more important for discontinuous fields (e.g. ring-and-plug) or regions of high strain gradient. Future implementations of this technique could implement slit width into the measurement model, treating each measurement as an area rather than line integral. 2. 2.

Intensity variation across the beam could be normalised in this case, but might be unavoidable at other neutron sources. Further improvements to the reconstruction algorithm could be made by implementing a weighting due to intensity in the measurement model. This weighting could either then be measured at the source in question or determined via hyperparameter optimisation. 3. 3.

Intensity variation due to sample texture was also not found to play a significant role in this experiment, but again could be implemented as an additional weighting in the measurement model and either measured or optimised for. 4. 4.

For the validation measurements, $d_{0}$ was calculated from the three measured strain components and the plane-stress condition. An average, constant $d_{0}$ was assumed for the tomography measurements. In some cases — particularly with unannealed/preprocessed samples, $d_{0}$ variation can result in significant pseudostrains. It is possible to estimate $d_{0}$ alongside the strains using the GP technique and such an extension to the reconstruction algorithm may be useful in some cases — see [16] for more information.

Extension of the algorithm presented here to 2D and 3D strain fields is conceptually straightforward, though practical issues may arise. A 2D implementation could involve taking projections of a sample by sweeping it past the secondary slit, then rotating the sample about it’s centre and repeating the process. A 2D reconstruction algorithm would likely resemble that in [11], with modifications made for attenuation and measurement direction.

While the results of this study found that diffraction tomography reconstructions can achieve results in agreement with conventional ‘pointwise’ measurements in equal or less beamtime, further work is required to show that this technique is viable for 2D systems, particularly with larger samples where attenuation of the beam may prove a limiting factor.

That being said, preliminary results suggest that — particularly for better scattering materials such as steel — projections of a similar confidence could be had in as little as 2 hours, compared to 7 hours to collect three components with conventional pointwise scanning (using a 1 mm slit width in both cases).

If a method for collecting entire projections at once is developed (e.g. using parallel collimators and a pixelated detector on a time-of-flight instrument), projections could be obtained in the time it presently takes to perform a single diffraction tomography measurement.

One of the more promising outcomes of this demonstration is the prospect of combined neutron transmission-diffraction tomography. In this case, both measurements could be performed concurrently and then processed by single reconstruction algorithm that takes into account their relative uncertainty.

Extension of this method to 2D systems, time-of-flight instruments and integration with existing transmission reconstruction algorithms forms a natural basis for future work.

5 Conclusion

In this paper, we proposed a new technique for strain measurement that we call neutron diffraction tomography. This technique draws on the recent activity in the field of Bragg-edge transmission strain tomography and aims to achieve similar results using conventional strain scanning geometry and instrumentation.

Reconstructions from simulated and experimental data were achieved and showed close agreement with validation measurements for axisymmetric systems with similar beamtime requirements.

We believe the method is viable. Future work involves extension to arbitrary 2D systems, where it may also be valuable in a combined transmission-diffraction strain reconstruction algorithm.

6 Acknowledgements

This work was supported by the Australian Research Council through a Discovery Project Grant No. DP17010 2324. Access to the KOWARI instrument was made possible through the ANSTO user-access program (Proposal No. 7318). The authors would also like to thank AINSE Limited for providing financial assistance (PGRA) and support to enable work on this project.

Appendix A Measurement Model Derivations

With respect to the coordinate system shown in Figure 3, we propose a measurement model of the form:

[TABLE]

where, $\boldsymbol{\hat{\kappa}}^{\mathsf{T}}\boldsymbol{\epsilon}(s)\boldsymbol{\hat{\kappa}}$ provides the average normal component of strain aligned with the measurement direction $\boldsymbol{\hat{\kappa}}$ in Cartesian coordinates by:

[TABLE]

We assume that the $y$ axis of our right-handed coordinate system is aligned with the path of diffracted neutron rays (in the direction of the detector). Thus, for a given diffraction angle $2\vartheta$, we have a measurement direction vector $\boldsymbol{\hat{\kappa}}$ given by:

[TABLE]

Converting strains from Cartesian to polar coordinates gives:

[TABLE]

In terms of the two nonzero in-plane components of strain in polar coordinates, we can then write:

[TABLE]

Where:

[TABLE]

and:

[TABLE]

The measurement can be separated into two integrals:

[TABLE]

And written as a linear transformation of the strain tensor:

[TABLE]

For future working, we define this transform $\boldsymbol{M}$:

[TABLE]