A Bayesian Approach to Triaxial Strain Tomography from High-energy X-ray Diffraction

TL;DR

This paper introduces a Bayesian Gaussian process model for 3D strain field tomography using high-energy X-ray diffraction, enabling rapid, high-resolution measurements with minimal rotation, and demonstrates its effectiveness through simulations.

Contribution

It presents a novel GP-based method for 3D strain tomography from high-energy X-ray data, incorporating static equilibrium constraints and achieving single-axis rotation reconstruction.

Findings

Effective in rejecting Gaussian noise in simulations

Capable of 3D strain reconstruction with a single rotation axis

Builds on previous 2D neutron measurement methods

Abstract

Diffraction of high-energy X-rays produced at synchrotron sources can provide rapid strain measurements, with high spatial resolution, and good penetrating power. With an uncollimated diffracted beam, through thickness averages of strain can be measured using this technique, which poses an associated rich tomography problem. This paper proposes a Gaussian process (GP) model for three-dimensional strain fields satisfying static equilibrium and an accompanying algorithm for tomographic reconstruction of strain fields from high-energy X-ray diffraction. We present numerical evidence that this method can achieve triaxial strain tomography in three-dimensions using only a single axis of rotation. The method builds upon recent work where the GP approach was used to reconstruct two-dimensional strain fields from neutron based measurements. A demonstration is provided from simulated data, showing the method is capable of rejecting realistic levels of Gaussian noise.

Equations

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

$$\begin{array}{rl}{I}_{ϵ}& =\frac{1}{L}\underset{0}{\overset{L}{\int }}{\hat{\mathbf{\kappa }}}^{\mathsf{T}}\mathbf{ϵ}(\mathbf{p}+\hat{\mathbf{n}}s)\hat{\mathbf{\kappa }}\mathrm{d}s+e,\end{array}$$

$$J_{ϵ}=\frac{1}{L}{\int }_0^L{\Pi }_{\hat{\mathbf{n}}}ϵ(\mathbf{p}+\hat{\mathbf{n}}s){\Pi }_{\hat{\mathbf{n}}},$$

$${\Pi }_{\hat{\mathbf{n}}}=\mathbf{I}-\hat{\mathbf{n}}{\hat{\mathbf{n}}}^{\mathsf{T}},$$

$${\Pi }_{\hat{\mathbf{n}}}ϵ{\Pi }_{\hat{\mathbf{n}}}=\left[\begin{array}{ccc}0& 0& 0\\ 0& {ϵ}_{yy}& {ϵ}_{yz}\\ 0& {ϵ}_{yz}& {ϵ}_{zz}\end{array}\right].$$

$$p(\mathbf{ϵ}|y)=\frac{p(y|\mathbf{ϵ})p(\mathbf{ϵ})}{p(y)}.$$

$$f(\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(m(\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime })\right),$$

$$\left[\begin{array}{c}f({\mathbf{x}}_1)\\ ⋮\\ f({\mathbf{x}}_N)\end{array}\right]\sim \mathcal{N}\left(\mathbf{\mu },\mathbf{K}\right)\text{where}\mu =\left[\begin{array}{c}m({\mathbf{x}}_1)\\ ⋮\\ m({\mathbf{x}}_N)\end{array}\right],$$

$$\mathbf{K}=\left[\begin{array}{ccc}k({\mathbf{x}}_1,{\mathbf{x}}_1)& \cdots & k({\mathbf{x}}_1,{\mathbf{x}}_N)\\ ⋮& \ddots & ⋮\\ k({\mathbf{x}}_N,{\mathbf{x}}_1)& \cdots & k({\mathbf{x}}_N,{\mathbf{x}}_N)\end{array}\right].$$

$$\mathbf{ϵ}(\mathbf{x})=\left[\begin{array}{ccc}{ϵ}_{xx}(\mathbf{x})& {ϵ}_{xy}(\mathbf{x})& {ϵ}_{xz}(\mathbf{x})\\ {ϵ}_{xy}(\mathbf{x})& {ϵ}_{yy}(\mathbf{x})& {ϵ}_{yz}(\mathbf{x})\\ {ϵ}_{xz}(\mathbf{x})& {ϵ}_{yz}(\mathbf{x})& {ϵ}_{zz}(\mathbf{x})\end{array}\right],$$

$$\overline{\mathbf{ϵ}}\sim \mathcal{G}\mathcal{P}\left(0,{\mathbf{K}}_{ϵ}(\mathbf{x},{\mathbf{x}}^{\prime })\right),$$

$$\mathbf{\sigma }(\mathbf{x})=\nabla \times \mathbf{\Phi }(\mathbf{x})\times \nabla ,$$

$$\mathbf{\sigma }(x,y,z)=\left[\begin{array}{ccc}{\sigma }_{xx}(\mathbf{x})& {\sigma }_{xy}(\mathbf{x})& {\sigma }_{xz}(\mathbf{x})\\ {\sigma }_{xy}(\mathbf{x})& {\sigma }_{yy}(\mathbf{x})& {\sigma }_{yz}(\mathbf{x})\\ {\sigma }_{xz}(\mathbf{x})& {\sigma }_{yz}(\mathbf{x})& {\sigma }_{zz}(\mathbf{x})\end{array}\right],$$

$$\mathbf{\Phi }(\mathbf{x})=\left[\begin{array}{ccc}{\Phi }_1(\mathbf{x})& {\Phi }_4(\mathbf{x})& {\Phi }_5(\mathbf{x})\\ {\Phi }_4(\mathbf{x})& {\Phi }_2(\mathbf{x})& {\Phi }_6(\mathbf{x})\\ {\Phi }_5(\mathbf{x})& {\Phi }_6(\mathbf{x})& {\Phi }_3(\mathbf{x})\end{array}\right].$$

$$\bar{\mathbf{\Phi }}=\left[\begin{array}{c}{\Phi }_1\\ {\Phi }_2\\ {\Phi }_3\\ {\Phi }_4\\ {\Phi }_5\\ {\Phi }_6\end{array}\right],\bar{\mathbf{\sigma }}=\left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{xy}\\ {\sigma }_{xz}\\ {\sigma }_{zz}\end{array}\right],\overline{\mathbf{ϵ}}=\left[\begin{array}{c}{ϵ}_{xx}\\ {ϵ}_{yy}\\ {ϵ}_{zz}\\ {ϵ}_{xy}\\ {ϵ}_{xz}\\ {ϵ}_{yz}\end{array}\right],$$

$$\left[\begin{array}{c}{ϵ}_{xx}\\ {ϵ}_{yy}\\ {ϵ}_{zz}\\ {ϵ}_{xy}\\ {ϵ}_{xz}\\ {ϵ}_{yz}\end{array}\right]=\frac{1}{E}\underset{\mathbf{H}}{\underbrace{\left[\begin{array}{cccccc}-\nu & 1& -\nu & 0& 0& 0\\ -\nu & -\nu & 1& 0& 0& 0\\ 0& 0& 0& 1+\nu & 0& 0\\ 0& 0& 0& 0& 1+\nu & 0\\ 0& 0& 0& 0& 0& 1+\nu \end{array}\right]}}\left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{xy}\\ {\sigma }_{xz}\\ {\sigma }_{zz}\end{array}\right],$$

$$\overline{\mathbf{ϵ}}=\mathbf{H}{\mathcal{B}}^{\mathbf{x}}\bar{\mathbf{\Phi }},$$

$${\mathcal{B}}^{\mathbf{x}}=\left[\begin{array}{cccccc}0& \frac{{\partial }^2}{\partial z^2}& \frac{{\partial }^2}{\partial y^2}& 0& 0& -2\frac{{\partial }^2}{\partial y\partial z}\\ \frac{{\partial }^2}{\partial z^2}& 0& \frac{{\partial }^2}{\partial x^2}& 0& -2\frac{{\partial }^2}{\partial x\partial z}& 0\\ \frac{{\partial }^2}{\partial y^2}& \frac{{\partial }^2}{\partial x^2}& 0& -2\frac{{\partial }^2}{\partial x\partial y}& 0& 0\\ 0& 0& -\frac{{\partial }^2}{\partial x\partial y}& -\frac{{\partial }^2}{\partial z^2}& \frac{{\partial }^2}{\partial y\partial z}& \frac{{\partial }^2}{\partial x\partial z}\\ -\frac{{\partial }^2}{\partial y\partial z}& 0& 0& \frac{{\partial }^2}{\partial x\partial z}& \frac{{\partial }^2}{\partial x\partial y}& -\frac{{\partial }^2}{\partial x^2}\\ 0& -\frac{{\partial }^2}{\partial x\partial z}& 0& \frac{{\partial }^2}{\partial y\partial z}& -\frac{{\partial }^2}{\partial y^2}& \frac{{\partial }^2}{\partial x\partial y},\end{array}\right]$$

$$\bar{\mathbf{\Phi }}(\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(\underset{{\mathbf{m}}_{\Phi }(\mathbf{x})}{\underbrace{\left[\begin{array}{c}{m}_1\\ m_2\\ m_3\\ m_4\\ m_5\\ m_6\end{array}\right]}},\underset{{\mathbf{K}}_{\Phi }(\mathbf{x},{\mathbf{x}}^{\prime })}{\underbrace{\left[\begin{array}{cccccc}{k}_1& 0& 0& 0& 0& 0\\ 0& k_2& 0& 0& 0& 0\\ 0& 0& k_3& 0& 0& 0\\ 0& 0& 0& k_4& 0& 0\\ 0& 0& 0& 0& k_5& 0\\ 0& 0& 0& 0& 0& k_6\end{array}\right]}}\right),$$

$$\begin{array}{rl}\overline{\mathbf{ϵ}}& \sim \mathcal{G}\mathcal{P}\left(\mathbf{H}{\mathcal{B}}^{\mathbf{x}}{\mathbf{m}}_{\Phi },\mathbf{H}{\mathcal{B}}^{\mathbf{x}}{\mathbf{K}}_{\Phi }(\mathbf{x},{\mathbf{x}}^{\prime }){\mathcal{B}}^{{\mathbf{x}}^{\prime }}{}^{\mathsf{T}}{\mathbf{H}}^{\mathsf{T}}\right)\\ & =\mathcal{G}\mathcal{P}\left(0,{\mathbf{K}}_{ϵ}(\mathbf{x},{\mathbf{x}}^{\prime })\right).\end{array}$$

$$\begin{array}{rl}{I}_{ϵ}(\mathbf{\eta })& =\frac{1}{L}\underset{0}{\overset{L}{\int }}\bar{\mathbf{\kappa }}\overline{\mathbf{ϵ}}(\mathbf{p}+\hat{\mathbf{n}}s)\mathrm{d}s+e\\ I_{ϵ}(\mathbf{\eta })& ={\mathcal{L}}^{\mathbf{x}}(\mathbf{\eta })\overline{\mathbf{ϵ}}(\mathbf{x})+e\end{array}$$

$$\left[\begin{array}{c}{\mathbf{I}}_{ϵ}\\ {\overline{\mathbf{ϵ}}}_{\ast }\end{array}\right]=\mathcal{N}\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{\mathbf{K}}_I+{\sigma }_n^2\mathbf{I}& {\mathbf{K}}_{\ast }^{\mathsf{T}}\\ {\mathbf{K}}_{\ast }& {\mathbf{K}}_{ϵ}\end{array}\right]\right),$$

$$\begin{array}{rl}({\mathbf{K}}_{\ast }{)}_i& ={\mathbf{K}}_{ϵ}{\mathcal{L}}^{{\mathbf{x}}^{\prime }}({\mathbf{\eta }}_i{)}^{\mathsf{T}}\\ & =\frac{1}{L_i}\underset{0}{\overset{L_i}{\int }}{\mathbf{K}}_{ϵ}({\mathbf{x}}_{\ast },{\mathbf{p}}_i+{\hat{\mathbf{n}}}_i{s}^{\prime }){\bar{\mathbf{\kappa }}}_i^{\mathsf{T}}\mathrm{d}{s}^{\prime }\end{array}$$

$$\begin{array}{rl}({\mathbf{K}}_I{)}_{ij}& ={\mathcal{L}}^{\mathbf{x}}({\mathbf{\eta }}_i){\mathbf{K}}_{ϵ}{\mathcal{L}}^{{\mathbf{x}}^{\prime }}({\mathbf{\eta }}_j{)}^{\mathsf{T}}\\ & =\frac{1}{L_i{L}_j}\underset{0}{\overset{L_i}{\int }}\underset{0}{\overset{L_j}{\int }}{\bar{\mathbf{\kappa }}}_i{\mathbf{K}}_{ϵ}({\mathbf{p}}_i+{\hat{\mathbf{n}}}_{i}s,{\mathbf{p}}_j+{\hat{\mathbf{n}}}_j{s}^{\prime }){\bar{\mathbf{\kappa }}}_j^{\mathsf{T}}\mathrm{d}{s}^{\prime }\mathrm{d}s\end{array}$$

$${\overline{\mathbf{ϵ}}}_{\ast }\sim \mathcal{N}({\mathbf{\mu }}_{\overline{ϵ}|{\mathbf{I}}_{ϵ}},{\mathbf{K}}_{{ϵ}_{\ast }|{\mathbf{I}}_{ϵ}})$$

$$\begin{array}{rl}{\mathbf{\mu }}_{{\overline{ϵ}}_{\ast }|{\mathbf{I}}_{ϵ}}& ={\mathbf{K}}_{\ast }({\mathbf{K}}_I+{\sigma }_n^2\mathbf{I}{)}^{-1}{\mathbf{I}}_{ϵ}\\ {\mathbf{K}}_{{ϵ}_{\ast }|{\mathbf{I}}_{ϵ}}& ={\mathbf{K}}_{ϵ}({\mathbf{x}}_{\ast },{\mathbf{x}}_{\ast })-{\mathbf{K}}_{\ast }({\mathbf{K}}_I+{\sigma }_n^2\mathbf{I}{)}^{-1}{\mathbf{K}}_{\ast }^{\mathsf{T}}\end{array}$$

$$\begin{array}{rl}{k}_i(\mathbf{x},{\mathbf{x}}^{\prime })& \approx \sum _{k=1}^m{\varphi }_{i,k}(\mathbf{x}){\Sigma }_{pi,kk}{\varphi }_{i,k}({\mathbf{x}}^{\prime })\\ & ={\mathbf{\varphi }}_i(\mathbf{x}){\Sigma }_{pi}{\mathbf{\varphi }}_i({\mathbf{x}}^{\prime }{)}^{\mathsf{T}},\end{array}$$

$$\begin{split}\phi_{i,k}(x)&=\frac{1}{L_{x}L_{y}L_{z}}\sin(\lambda_{xk}(x+L_{x}))\\ &\hskip 28.45274pt\sin(\lambda_{yk}(y+L_{y}))\sin(\lambda_{zk}(z+L_{z})),\\ \Sigma_{pi,k}&=\int K(\mathbf{r})\exp({\color[rgb]{0,0,0}\mathrm{-i}}\boldsymbol{\lambda}^{\mathsf{T}}\mathbf{r})\,\mathrm{d}r.\end{split}$$

$${\Sigma }_{pi,kk}={\sigma }_f^2(2\pi {)}^{\frac{3}{2}}{l}_x{l}_y{l}_z\exp \left(-\frac{1}{2}\left(l_x^2{\lambda }_{xk}^2+l_y^2{\lambda }_{yk}^2+l_z^2{\lambda }_{zk}^2\right)\right),$$

$$\begin{array}{rl}{\mathbf{K}}_{\Phi }& \approx {\mathbf{\varphi }}_{\Psi }{\mathbf{\Sigma }}_{\Psi }{\mathbf{\varphi }}_{\Psi }^{\mathsf{T}},\\ {\mathbf{\varphi }}_{\Psi }& ={\left[\begin{array}{c}{\mathbf{\varphi }}_1\\ {\mathbf{\varphi }}_2\\ {\mathbf{\varphi }}_3\\ {\mathbf{\varphi }}_4\\ {\mathbf{\varphi }}_5\\ {\mathbf{\varphi }}_6\end{array}\right]}^{\mathsf{T}},{\mathbf{\Sigma }}_{\Psi }=\left[\begin{array}{cccccc}{\Sigma }_{p1}& 0& 0& 0& 0& 0\\ 0& {\Sigma }_{p2}& 0& 0& 0& 0\\ 0& 0& {\Sigma }_{p3}& 0& 0& 0\\ 0& 0& 0& {\Sigma }_{p4}& 0& 0\\ 0& 0& 0& 0& {\Sigma }_{p5}& 0\\ 0& 0& 0& 0& 0& {\Sigma }_{p6}\end{array}\right]\end{array}$$

$$\begin{array}{rl}{\mathbf{K}}_{\mathbf{ϵ}}& \approx {\mathbf{\varphi }}_{ϵ}({\mathbf{x}}_{\ast }){\mathbf{\Sigma }}_{\Psi }{\mathbf{\varphi }}_{ϵ}({\mathbf{x}}_{\ast }{)}^{\mathsf{T}}\\ ({\mathbf{K}}_{\ast }{)}_i& \approx {\mathbf{\varphi }}_{ϵ}({\mathbf{x}}_{\ast }){\mathbf{\Sigma }}_{\Psi }{\mathbf{\varphi }}_{\mathbf{I}}({\mathbf{\eta }}_i{)}^{\mathsf{T}}\\ ({\mathbf{K}}_{\mathbf{I}}{)}_{ij}& \approx {\mathbf{\varphi }}_{\mathbf{I}}({\mathbf{\eta }}_i){\Sigma }_{\Psi }{\mathbf{\varphi }}_{\mathbf{I}}({\mathbf{\eta }}_j{)}^{\mathsf{T}}\end{array}$$

Full text

††journal: ArXiv

1 Introduction

Diffraction of X-rays and neutrons allow the study of mechanical stress and strain within crystalline solids [10, 15, 6, 5]. These techniques revolve around Bragg’s law; $\lambda=2d\sin\theta$, whereby variations in the lattice spacing $d$, due to elastic strain can be observed by changes in the scattering angle, $\theta$, of the diffracted radiation with wavelength $\lambda$. These variations are related to the average elastic strain within the scattering volume according to

[TABLE]

where $d_{0}$ is the lattice spacing in a strain free sample. Both neutron and conventional X-ray diffraction have limitations; lab-based X-rays sample only a shallow surface layer (typically a few microns), while the characteristic low intensities of neutrons gives rise to long acquisition times and spatial resolutions of $1\text{mm}$ or larger.

Modern X-ray synchrotron sources can provide very intense narrow beams of highly penetrating X-ray photons [22]. These high-energy X-rays can provide strain measurements with beam spot sizes as small as $1\mu\text{m}$ and can have path lengths of many centimetres, even in steel. Although high-energy X-ray diffraction can be used to study dynamic strains (e.g. [14]), in this paper the problem is restricted to steady state. The restriction to steady-state allows the static equilibrium constraints to be assumed as prior information.

A particular feature of high energy X-rays is that the scattering angle is typically small $2\theta<15^{\circ}$ [22]. This has two implications for the study of strain fields: firstly, that the normal component of strain measured is almost perpendicular to the incident beam; at an angle of $90^{\circ}-\theta$; secondly, the gauge volume is elongated in the direction of the incident beam (typically by a factor of 10 or more).

The ability to select small beam cross-sections allows strain profiles [3] or two dimensional strain maps [22] to be readily obtained in the plane perpendicular to the incident beam (see Figure 1 for the measurement geometry). However, the resolution along the incident beam is degraded by a factor of 10 or more due to the elongation of the scattering volume. This does not always present a problem, for example in [3] the through thickness strain variation in the direction of the incident beam was known to be small, and averaging in this direction provided a good measure of of strain.

To overcome these limitations, a different approach is presented in [11] and [12] where determining the strain field is considered a rich tomography problem from a series of through thickness measurements. In this setting, the goal is to reconstruct the higher order (two- or three-dimensional) distribution of unknown strain from a set of lower dimension (one- or two-dimensional) projections. In [11], the axisymmetric strain within a quenched cylinder was reconstructed, providing an initial demonstration of this conceptual approach. This approach was extended in [12], where reconstruction of the longitudinal strain within a zirconia dental prostheses was achieved by posing the problem in a form suitable for conventional computed tomography algorithms. In this experiment, the sample was rotated about a single axis while through-thickness averages of the out-of-plane normal strain were recorded. The resulting scalar tomography problem was then solved using conventional back-projection techniques.

This problem is also studied by Lionheart and Withers[13], where the transverse ray transform [18] is given as a possible model for the high-energy X-ray strain measurements [13]. This assumes that the direction of measured strain is exactly perpendicular to the beam, i.e. that the angle of diffraction is close to zero. Further, It is proposed that full triaxial strain field tomography can be done as a series of ‘regular’ tomography problems. Here, each regular tomography involves rotating the sample about a single axis and reconstructing the component of strain in the direction of the axis of rotation. By performing six of these experiments about different axes of rotation the full triaxial strain field could be determined. The conditions on the the axes chosen to allow for the the triaxial strain to be reconstructed are given.

This work is extended in Desai and Lionheart[4] through the presentation of an explicit inversion formula for the transverse ray transform. In general, it was shown that this formula allows for the reconstruction of strain from high-energy X-ray measurements made around three axes of rotation. It was also shown that if compatibility can be assumed then only two axes of rotation are required. Decreasing the number of axes that the sample has to be rotated about reduces the complexity of the experiment set up and the time required to run the experiment.

In this paper, a Gaussian process based approach to reconstruct the triaxial strain field from these measurements is presented. This approach enforces that the reconstruction to satisfies the static equilibrium constraints, and is applied to the problem of reconstructing the strain field using single axis tomography with promising results in simulation. Recent work in related fields has demonstrated the use of Gaussian processes (GPs) to model and reconstruct two-dimensional strain fields from time-of-flight neutron transmission measurements [9, 7] and neutron diffraction measurements [8].

1.1 Contribution

This paper makes the following contributions:

A Gaussian process model for steady-state, triaxial elastic strain fields in three-dimensions. This is a non-trivial extension of the model presented in Jidling et al [9] which was restricted to biaxial strain fields in two dimensions under the assumption of plane stress or plane strain. 2. 2.

A method for applying this model to reconstruct the strain field from high-energy X-ray strain measurements. Results indicate that by assuming equilibrium this approach can reconstruct a triaxial strain field using a single axis of rotation. This could not be achieved using existing methods.

2 High Energy X-ray Strain Measurement

This section provides a brief summary of strain measurement using high energy X-ray diffraction. This process is not the primary focus of this work, however the summary provides details pertinent to model the relationship between the measurements and the strain field. A more detailed description can be found in [22] and details of geometric corrections required due to rotating the sample causing changes in the sample to detector distance are describe in [12]. The summary is as follows.

The incident beam with direction $\hat{\mathbf{n}}$ is diffracted within the sample according to Bragg’s law at an angle of $2\theta$ forming a shallow cone. The intensity of the diffracted beam is recorded at a detector, with the peak intensity’s forming a Debye-Scherrer ring, see Figure 1. Due to the small diffraction angle, the distance to the detector, $D$, is much greater than nominal sample dimensions.

For a polychromatic X-ray beam, a diffraction pattern can be fit to a segment of the Debye-Scherrer ring providing a measurement of the average normal strain in the direction$\kappa$ of the form (1); with a $10^{\circ}$ being suitable for strain measurements [12]. Two possible segments and are shown in Figure 1 with corresponding measurement directions $\hat{\boldsymbol{\kappa}}_{1}$ and $\hat{\boldsymbol{\kappa}}_{2}$, respectively. Multiple segments can be analysed from each Debye-Scherrer ring to give measurements of the normal strain in different directions. Although using the relative shift of the diffraction pattern has been done in practice, Lionheart and Withers[13] shows that, theoretically, a particular moment of the diffraction pattern should be used instead.

In this work, we consider the diffracted beam to be left uncollimated—which aligns with the work in [12], and the incident beam to be collimated to give a spot size of $h\times h$. This measurement geometry is shown in Figure 2 and gives the scattering volume’s length as the path length through the sample, $L$.

A measurement model can be formulated using the following reasoning. The measurement corresponds to an average of normal strain in direction $\hat{\boldsymbol{\kappa}}$ over the gauge volume. This average can be represented by a volume integral divided by the volume. However, since spot sizes as small as $1\mu\text{m}$ are achievable and that typical path lengths can be several orders of magnitude larger, it is reasonable to assume that there is no variation in the perpendicular to the direction of the beam. Under this assumption, the volume integral can be reduced to a line integral in the beam direction;

[TABLE]

where $e\sim\mathcal{N}(0,\sigma_{n}^{2})$, $\boldsymbol{\epsilon}(\mathbf{x})$ is the strain tensor field inside the sample and $\mathbf{p}$ is the initial intersection between the beam and the sample.

By analysing multiple segments from each Debye-Scherrer ring, each ray provides information about multiple components of strain. For example, assuming that the measurement direction is exactly perpendicular to the beam direction. Then, choosing at least three unique measurement directions, each ray provides information about the strain projected onto the plane perpendicular to the beam. That is, a beam with direction $\hat{\mathbf{n}}=[1,0,0]^{\mathsf{T}}$ predominantly provides information about $\epsilon_{yy}$, $\epsilon_{yz}$, and $\epsilon_{zz}$.

If we make the additional assumption that the measurement directions are perpendicular to the incident beam, then this line of reasoning also leads to the transverse ray transform (TRT) which was proposed by [13] to model high-energy X-ray strain measurements. With the addition of noise, the TRT model for the measurement is given by

[TABLE]

where

[TABLE]

and $\mathbf{I}$ is the identity matrix. Here, $\Pi_{\hat{\mathbf{n}}}\epsilon\Pi_{\hat{\mathbf{n}}}$ is the projection of the strain field onto the plane perpendicular to $\hat{\mathbf{n}}$. For example, consider that $\hat{\mathbf{n}}=[1,0,0]^{\mathsf{T}}$, then

[TABLE]

Although the TRT represents a more compact measurement model, it requires assuming the measurement directions are perpendicular to the beam direction, which is not quite the case in practice. Hence, (2) has been used as the measurement model in this work.

3 3D Strain field Reconstruction

In this work, the strain field is reconstructed using the framework of Bayesian inference. The strain field is modelled as having a probability distribution, this does not mean that the strain field is random, instead it represents our uncertainty in its values. A prior probability distribution is assigned to the strain field, $p(\boldsymbol{\epsilon})$, that represents any knowledge we have before the inclusion of measurements, $y$. This distribution is then updated by the inclusion of measurements, which are assigned a likelihood $p(y|\boldsymbol{\epsilon})$, using Bayes’ rule to give a posterior distribution of the strain:

[TABLE]

As the strain field must satisfy equilibrium, the prior distribution should incorporate this knowledge. In two dimensions this can be achieved by modelling the Airy’s stress function by a Gaussian process (GP) from which the two-dimensional strain field can be defined under an assumption of plane stress or plane strain [9]. This method was successfully used to reconstruct two-dimensional strain fields from neutron transmission measurements [9, 7] and neutron diffraction measurements [8].

The following sections present a generalisation of this procedure to three-dimensions and application to reconstructing the strain field from from high energy X-ray diffraction measurements. A brief overview of GPs is provided in Section 3.1. A Gaussian process prior suitable for modelling three-dimensional strain fields is designed in Section 3.2. A likelihood model for high energy X-ray strain measurements is then defined in Section 3.3, allowing the strain field to be reconstructed. For this work, the material should have a randomly distributed polycrystalline structure with no preferred crystal orientation (i.e. no texture) such that the bulk material is isotropic.

3.1 Gaussian Process

This section gives a brief description of GPs, for a more detailed description see Rasmussen and Williams[16]. A GP is a stochastic process suitable for modelling spatially correlated functions and can be viewed as a distribution over functions;

[TABLE]

where $m(\mathbf{x})=\mathbb{E}\left[f(\mathbf{x})\right]$ is the mean function and $k(\mathbf{x},\mathbf{x}^{\prime})=\mathbb{E}\left[\left(m(\mathbf{x})-f(\mathbf{x})\right)\left(m(\mathbf{x})-f(\mathbf{x})\right)^{\mathsf{T}}\right]$ is the covariance function. The choice of covariance function governs the characteristics of the functions in this distribution—such as their smoothness. Many choices exist for the covariance function and a good summary is available elsewhere [16].

A Gaussian process is a generalisation of a multivariate Gaussian in the sense that function values sampled at a finite number of inputs $\mathbf{x}_{1},\dots,\mathbf{x}_{N}$ are Gaussian distributed;

[TABLE]

If a GP prior is used and the measurement likelihood is Gaussian then posterior distribution given by Bayes’ rule has a closed form and can be computed using standard Gaussian conditioning.

3.2 Three-dimensional Strain Field GP

Here, we design a GP prior for strain in three-dimensions that intrinsically satisfies static equilibrium. The strain field is defined as the symmetric tensor

[TABLE]

where $\mathbf{x}=[x\ \ y\ \ z]^{\mathsf{T}}$ are the spatial coordinates.

This GP will be of the form

[TABLE]

where $\bar{\boldsymbol{\epsilon}}$ is a vector of the $6$ unique components of the strain field. Here, the covariance function for the strain field, $\mathbf{K}_{\epsilon}(\mathbf{x},\mathbf{x}^{\prime})$ will be designed to ensure that all strain fields belonging to this GP satisfy equilibrium. This is done by specifying a GP prior for a set of potentials known as the Beltrami stress functions. Having done this, a GP prior for the strain field can be derived using the equilibrium equations and Hooke’s law. The details of this derivation are as follows.

The Beltrami stress functions [1] allow a complete solution to the equilibrium equations in three-dimensions to be written as [20, 17, 2];

[TABLE]

where $\boldsymbol{\nabla}=\begin{bmatrix}\frac{\partial}{\partial x}&\frac{\partial}{\partial y}&\frac{\partial}{\partial z}\end{bmatrix}^{\mathsf{T}}$ and $\boldsymbol{\sigma}(x,y,z)$ is the symmetric stress tensor;

[TABLE]

and $\boldsymbol{\Phi}(x,y,z)$ is the Beltrami stress tensor consisting of six unique scalar potential fields;

[TABLE]

To improve readability we introduce the following vectorised notation;

[TABLE]

where the spatial coordinates $\mathbf{x}$, are omitted for brevity. Using this notation we can write Hooke’s law which relates the stress and strain fields as

[TABLE]

where $\mathbf{H}$ is the compliance matrix for isotropic materials, $\nu$ is Poisson’s ratio, and $E$ is Young’s modulus. The scaling $\frac{1}{E}$ can be safely neglected in the methods implementation to provide better numerical scaling. Finally, we can write the mapping from the Beltrami stress functions to the strain field in this vectorised form as

[TABLE]

where

[TABLE]

and the superscript is used to denote the set of spatial coordinates on which the operator acts.

As this mapping is linear it can be used to define a GP prior on the strain function [21, 16, 9]. Each component of the Beltrami tensor field, $\Phi_{i}(\mathbf{x})$, is assigned a Gaussian process prior with its own covariance function $k_{i}(\mathbf{x},\mathbf{x}^{\prime})$ — the squared-exponential is used in Section 5. The covariance function encodes spatial correlation which assumes the resulting strain field will have a degree of smoothness. Each component of the Beltrami tensor is considered independent as their is no prior information to suggest otherwise, and including correlation between the components would result in the strain field being restricted to meet an additional conditions beyond the desired equilibrium condition. This gives a GP for $\bar{\boldsymbol{\Phi}}(\mathbf{x})$ as

[TABLE]

where the shorthand $k_{i}=k_{i}(\mathbf{x},\mathbf{x}^{\prime})$ and $m_{i}=m_{i}(\mathbf{x})$ has been used.

The mapping (16) is applied to give a GP prior for the strain field that will ensure that any estimated strain field satisfies equilibrium;

[TABLE]

The covariance function for the strain field, $\mathbf{K}_{\epsilon}$, has correlations between the individual components of strain that ensure any estimated strain field will satisfy equilibrium. Here, without loss of generality [16], we chose the prior mean functions as zero.

3.3 Reconstruction from X-ray Strain Measurements

Here, we define the likelihood of the X-ray strain measurements $p(I_{\epsilon}|\boldsymbol{\epsilon})$. In vector form we can write the measurement model as

[TABLE]

where $\bar{\boldsymbol{\kappa}}=\begin{bmatrix}\kappa_{x}^{2}&\kappa_{y}^{2}&\kappa_{z}^{2}&2\kappa_{x}\kappa_{y}&2\kappa_{x}\kappa_{z}&2\kappa_{y}\kappa_{z}\end{bmatrix}$, $\boldsymbol{\eta}=\left\{\boldsymbol{\kappa},\hat{\mathbf{n}},\mathbf{p},L\right\}$, $e\sim\mathcal{N}(0,\sigma_{n}^{2})$, and $\mathcal{L}^{\mathbf{x}}(\boldsymbol{\eta})$ is a considered an operator that maps from the strain function $\bar{\boldsymbol{\epsilon}}(\mathbf{x})$ into the measurements $I_{\epsilon}(\boldsymbol{\eta})$. As this operator is linear, the joint distribution of the strain field at user specified location of interest, $\bar{\boldsymbol{\epsilon}}_{*}=\bar{\boldsymbol{\epsilon}}(\mathbf{x}_{*})$, and the measurements $\mathbf{I}_{\epsilon}=\begin{bmatrix}I_{\epsilon}(\boldsymbol{\eta}_{1})&\cdots&I_{\epsilon}(\boldsymbol{\eta}_{n})\end{bmatrix}^{\mathsf{T}}$ is Gaussian [21, 16, 9];

[TABLE]

where $\mathbf{K}_{*}$ is the cross covariance between the strains and the measurements and $\mathbf{K}_{I}$ is the covariance of the measurements.

The cross covariance between the strain $\bar{\boldsymbol{\epsilon}}_{*}$ and a measurement $I_{\epsilon}(\boldsymbol{\eta}_{i})$ is given by a single application of (20) to the strain fields covariance function;

[TABLE]

and the covariance between each pair of measurements, $I_{\epsilon}(\boldsymbol{\eta}_{i})$ and $I_{\epsilon}(\boldsymbol{\eta}_{j})$ is similarly given by

[TABLE]

The posterior distribution of the strain $\bar{\boldsymbol{\epsilon}}$ conditioned on the measurements $\mathbf{I}_{\epsilon}$ is

[TABLE]

where

[TABLE]

An analytic solution to the double integral is not known. However, computationally expensive numerical integration can be avoided by using an approximation scheme such as the one described in Section 4.1. It is worth noting that the extension to non-convex geometry is straight forward [9].

4 Implementation

Two practical aspects need to be considered when implementing the reconstruction algorithm; computational complexity, and hyperparameter selection.

4.1 Reducing Computation Complexity

The computational complexity of solving Equation (25) is twofold: firstly, the construction of $\mathbf{K}_{II}$ requires the evaluation of a double integral for every unique pair of measurements; secondly, the time required to invert $\mathbf{K}_{II}+\sigma_{n}^{2}\mathbf{I}$ scales with $\mathcal{O}(N^{3})$. X-ray strain tomography problems of the type considered here have a large number of measurements, motivating the use of an approximation scheme. Here, we consider the approximation scheme proposed in [19] which has previously been used for strain estimation [9, 8]. Using this scheme the covariance functions assigned to the Beltrami stress functions, $\mathbf{K}_{\Phi_{i}}$, can be approximated by a finite series of $m$ basis functions;

[TABLE]

where each $\phi_{i,k}(\mathbf{x})$ is a basis function and $\Sigma_{pi,kk}$ is its spectral density;

[TABLE]

Here, $\mathbf{r}=\mathbf{x}-\mathbf{x}^{\prime}$, and $\boldsymbol{\lambda}=[\lambda_{xj}\ \ \lambda_{yj}\ \ \lambda_{zj}]^{\mathsf{T}}$. For the squared-exponential covariance function used in Section 5 the spectral density is

[TABLE]

where $\sigma_{f}$, $l_{x}$, $l_{y}$, and $l_{z}$ are commonly referred to as hyper-parameters and are discussed in Section 4.2. The parameters $L_{x}$, $L_{y}$, $L_{z}$, and $\boldsymbol{\lambda}$ are analogous to the frequency and phase of the basis functions.

In this work, they were chosen so that the basis functions spanned a region where their spectral densities were greater than a predefined minimum threshold. This helps to ensure numerical stability while capturing the dominant modes of the reconstruction.

By concatenating the basis functions we can concisely express the approximation for the covariance function $\mathbf{K}_{\Phi}$;

[TABLE]

It is now a straight forward application of the mappings (16) and (20) to approximate the covariances required in Section 3 to reconstruct the strain field;

[TABLE]

This simplifies the problem to the calculation of $\boldsymbol{\phi}_{\epsilon}(\mathbf{x}_{*})$ and $\boldsymbol{\phi}_{I,i}(\boldsymbol{\eta}_{i})$, which are intuitively the basis functions for the strain field and the measurements, respectively, and only require a single application of the mappings. Complete expressions for these basis functions are found in Appendix A.

The posterior mean and covariance are approximated by

[TABLE]

This avoids forming the covariance matrices, reduces the complexity of the regression to $\mathcal{O}(Nm^{2})$ and also removes the need for numerical derivatives.

4.2 Hyperparameter Optimisation

The covariance functions, assigned to the potential functions, are characterised by their hyperparameters $\boldsymbol{\theta}$. For example, the squared-exponential covariance function has hyperparameters $\boldsymbol{\theta}=\{\sigma_{f},l_{x},l_{y},l_{z}\}$; where $\sigma_{f}$ encodes our prior uncertainty and the length scales $l_{x}$, $l_{y}$, and $l_{z}$ provide an assumption of smoothness. The hyperparmaters are selected by maximising the marginal log likelihood, $\log p(\mathbf{I}_{\epsilon}|\{\boldsymbol{\eta}_{i}\},\boldsymbol{\theta})$;

[TABLE]

where $\mathbf{K}_{II}$ is a function of $\boldsymbol{\theta}$. For the approximation in Section 4.1 $\mathbf{K}_{II}$ is replaced by $\boldsymbol{\phi}_{I}\Sigma_{p}\boldsymbol{\phi}_{I}^{\mathsf{T}}$, and expressions for the approximate log likelihood and its derivatives can be found in [8] or [9]. The parameters can thereafter be optimised using a gradient-based method, such as the BFGS algorithm in [23].

5 Simulation Demonstration

The method presented is demonstrated on a simulated 3D cantilevered beam example with the strain field given by a superposition of the Saint-Venant approximation to the strain field in the $xy$- and $yz$-planes. This is chosen as an appropriate example in the absence of a suitable experimental data set. The strain field is given by

[TABLE]

where geometry and loading is defined in Figure 3.

Measurements were simulated through these strain fields using Equation (20). The measurement geometry corresponded to rotating the sample about $z$ and for each angle using a $40\times 40$ grid of incident beams. A diffraction angle of $2\theta=10^{\circ}$ was used for the simulation, giving the angle between the incident beam and the direction of measured strain as $\alpha=85^{\circ}$. For each incident beam 36 strain directions are measured, corresponding to using $10^{\circ}$ segments from the Debye-Scherrer rings. The measurements were corrupted by zero-mean Gaussian noise with standard deviation $\sigma_{m}=1\times 10^{-4}$.

A convergence study as the number of rotation angles is increased was run and the results are shown in Figure 4. For each simulation the $n_{\psi}$ angles were chosen to be linearly spaced over $\frac{n_{\psi}}{n_{\psi}+1}180^{\circ}$. From equation 20 it is clear that two beams with travelling along the same path but with opposite directions provide measurements of almost identical components of the strain field. Therefore, it is not necessary to choose angles spanning $360^{\circ}$. The mean absolute error and the mean marginal standard deviation of the predicted strains are reported. The strong correlation between the mean marginal standard deviation and the mean absolute error suggests that the posterior covariance $\mathbf{K}_{\epsilon_{*}|\mathbf{I}}$ could inform the user as to the expected error in the reconstruction.

These results indicate that at least three rotation angles are required for an accurate reconstruction. This is true regardless of the number of incident beams used per rotation angle, and intuitively can be explained as each rotation angle the incident beams predominantly provide information about the strains lying in the plane perpendicular to the beams. Hence, three unique sets of these ‘in-plane’ strains need to be observed to recover the three-dimensional strain field.

The results of reconstructing from the measurement sets generated using three and twenty rotation angles are shown in Figure 5. The results show the reconstructed effective and hydrostatic strains along two cutting planes, which give a good indication of the overall accuracy of the reconstruction. The effective111The effective strain is the strain equivalent of von Mises stress sometimes referred to as the Mises effective strain., $\epsilon_{\text{eff}}$, and hydrostatic, $\epsilon_{\text{hyd}}$, strains are defined as

[TABLE]

Numerical simulations were also used to investigate the importance of $\alpha$ on the ability of the method to accurately reconstruct the the strain field. Comparing the accuracy when the strain field was reconstructed from 10 projections with $\alpha$ ranging from 90 to 85 degrees it was found that the specific angle was not important for the proposed method to be able to accurately reconstruct the strain field. For example, using $\alpha=90^{\circ}$ gave a mean relative error of $1.04\%$ compared to $0.99\%$ for $\alpha=85^{\circ}$. This supports the argument that it is the assumption of equilibrium rather than measurement direction that are not quite perpendicular to the incident beam that allows the proposed method to successfully reconstruct the simulated strain field using single-axis tomography.

In addition to the mean value for the reconstructed strain field, Equation 25 also provides a measure of the uncertainty. This can be used to determine where the reconstruction is least certain. For example, the marginal standard deviation of the reconstructed $\epsilon_{xx}$ component from $20$ projections is shown in Figure 6 for the cutting plane at $z=0\text{mm}$. The standard deviation is highest near the sample boundary indicating that this is where the largest errors are likely to be. Standard deviation maps of the other strain components are similar.

6 Conclusion

A Bayesian approach to triaxial strain tomography from high-energy X-ray diffraction measurements has been presented. This approach models the strain field as a Gaussian process such that the resulting reconstruction always satisfies equilibrium. From simulated measurements, this method was found to be capable of accurately reconstructing a full strain field in the presence of realistic levels of measurement noise. These numerical results suggest that it is possible to tomographically reconstruct a full triaxial strain distribution using single axis tomography by assuming equilibrium. While the squared-exponential covariance function proved suitable in this demonstration, other covariance functions, such as the Matérn or rational quadratic, may be more suitable for a particular strain field.

The next stage of this work will involve the planning and execution of an experimental demonstration of the technique. The ability to test the algorithm on simulated strain fields will provide guidance in this process. For example, for a given strain field the trade-off between measurement noise, number of projection angles, and spacing of measurements for each projection can be investigated.

7 Acknowledgements

This work is supported by the Australian Research Council through a Discovery Project Grant (DP170102324).

Appendix A Approximation Basis Functions

This section provides expressions for the basis functions given in (30) required for the approximation scheme used in Section 4.1.

As stated earlier, a set of bases for the strain field is defined by a linear mapping from a set of bases for the Beltrami stress functions;

[TABLE]

This is a linear combination of the partial derivatives of each Beltrami stress, i, function basis, k;

[TABLE]

where

[TABLE]

Applying the measurement mapping, $\mathcal{L}(\boldsymbol{eta}_{i})$ to the basis functions for the strain field gives a set of basis functions for the $i^{\text{th}}$ measurement;

[TABLE]

Which is a linear combination of the line integrals of the partial derivatives of each basis function. Defining

[TABLE]

the necessary components can be written as

[TABLE]