Neutron Transmission Strain Tomography for Non-Constant Stress-Free Lattice Spacing

TL;DR

This paper introduces a new neutron transmission strain tomography method that jointly reconstructs strain and stress-free lattice spacing fields, overcoming previous limitations of assuming a constant $d_0$, and demonstrates its effectiveness on simulated data.

Contribution

It presents a novel non-linear reconstruction algorithm for simultaneous strain and $d_0$ field estimation, incorporating equilibrium and boundary conditions.

Findings

Successful joint reconstruction of strain and $d_0$ fields on simulated data.

Method handles realistic noise levels effectively.

Ensures reconstructed strain fields satisfy physical equilibrium constraints.

Abstract

Recently, several algorithms for strain tomography from energy-resolved neutron transmission measurements have been proposed. These methods assume that the stress-free lattice spacing $d_0$ is a known constant limiting their application to the study of stresses generated by manufacturing and loading methods that do not alter this parameter. In this paper, we consider the more general problem of jointly reconstructing the strain and $d_0$ fields. A method for solving this inherently non-linear problem is presented that ensures the estimated strain field satisfies equilibrium and can include knowledge of boundary conditions. This method is tested on a simulated data set with realistic noise levels, demonstrating that it is possible to jointly reconstruct $d_0$ and the strain field.

Figures

Equations

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

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

$$y(\mathbf{\eta })=d(\mathbf{\eta })+e=\frac{1}{L}\underset{0}{\overset{L}{\int }}\bar{\mathbf{n}}\overline{\mathbf{ϵ}}(\mathbf{p}+\hat{\mathbf{n}}s)d_0(\mathbf{p}+\hat{\mathbf{n}}s)+d_0(\mathbf{p}+\hat{\mathbf{n}}s)\mathrm{d}s+e.$$

$$\begin{array}{rl}\frac{\partial }{\partial x}({ϵ}_{xx}+\nu {ϵ}_{yy})+\frac{\partial }{\partial y}(1-\nu ){ϵ}_{xy}& =0,\\ \frac{\partial }{\partial x}({ϵ}_{yy}+\nu {ϵ}_{xx})+\frac{\partial }{\partial y}(1-\nu ){ϵ}_{xy}& =0,\end{array}$$

$$0=\left[\begin{array}{ccc}{n}_{\perp 1}& n_{\perp 2}& 0\\ 0& n_{\perp 1}& n_{\perp 2}\end{array}\right]\left[\begin{array}{ccc}1& 0& -\nu \\ 0& 1+\nu & 0\\ -\nu & 0& 1\end{array}\right]\overline{\mathbf{ϵ}}({\mathbf{x}}_b),$$

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

$$\overline{\mathbf{ϵ}}(\mathbf{x})={\mathcal{V}}^{\mathbf{x}}\phi (\mathbf{x}),{\mathcal{V}}^{\mathbf{x}}=\left[\begin{array}{c}\frac{{\partial }^2}{\partial y^2}-\nu \frac{{\partial }^2}{\partial x^2}\\ -(1+\nu )\frac{{\partial }^2}{\partial x\partial y}\\ \frac{{\partial }^2}{\partial x^2}-\nu \frac{{\partial }^2}{\partial y^2}\end{array}\right],$$

$$\overline{\mathbf{ϵ}}(\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(0,{\mathcal{V}}^{\mathbf{x}}{k}_{\phi }(\mathbf{x},{\mathbf{x}}^{\prime }){\mathcal{V}}^{{\mathbf{x}}^{\prime }}{}^{\mathsf{T}}\right),$$

$$d_0(\mathbf{x})\sim \mathcal{G}\mathcal{P}({\bar{d}}_0,k_{d_0}(\mathbf{x},{\mathbf{x}}^{\prime })).$$

$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\sum _{j=1}^m{\varphi }_i(\mathbf{x})S({\mathbf{\lambda }}_j){\varphi }_j({\mathbf{x}}^{\prime }),$$

$$S(\mathbf{\omega })=\int k(\mathbf{r})e^{-i{\mathbf{\omega }}^{\mathsf{T}}\mathbf{r}}\mathrm{d}\mathbf{r},{\varphi }_j=\frac{1}{\sqrt{L_x{L}_y}}\sin ({\lambda }_{x,j}(x+L_x))\sin ({\lambda }_{y,j}(y+L_y)),$$

$$f(\mathbf{x})=\sum _{j=1}^m{\varphi }_j(\mathbf{x})w_j=\mathbf{\varphi }(\mathbf{x})\mathbf{w},w_j\sim \mathcal{N}({\mu }_j,S({\mathbf{\lambda }}_j)),$$

$$\begin{array}{rl}& {\overline{\mathbf{ϵ}}}_{\ast }(\mathbf{x})={\mathbf{\varphi }}_{ϵ}{\mathbf{w}}_{\phi },{\varphi }_{ϵ,j}(\mathbf{x})={\mathcal{V}}^{\mathbf{x}}{\varphi }_{\phi ,j}(\mathbf{x}),{\varphi }_{\phi ,j}=\frac{1}{\sqrt{L_{\phi x}{L}_{\phi y}}}\sin ({\lambda }_{\phi x,j}(x+L_{\phi x}))\sin ({\lambda }_{\phi y,j}(y+L_{\phi y})),\\ & d_{0\ast }(\mathbf{x})={\mathbf{\varphi }}_{d_0}{\mathbf{w}}_{d_0},{\varphi }_{d_0,k}(\mathbf{x})=\frac{1}{\sqrt{L_{d_{0}x}{L}_{d_{0}y}}}\sin ({\lambda }_{d_{0}x,k}(x+L_{d_{0}x}))\sin ({\lambda }_{d_{0}y,k}(y+L_{d_{0}y})),\end{array}$$

$$\begin{array}{rl}{y}_{\ast }& =\frac{1}{L}\underset{0}{\overset{L}{\int }}\bar{\mathbf{n}}\left(\sum _j\sum _k{\varphi }_{ϵ,j}(\mathbf{p}+\hat{\mathbf{n}}s)w_{\phi ,j}{\varphi }_{d_0,k}(\mathbf{p}+\hat{n}s)w_{d_0,k}\right)+\left(\sum _k{\varphi }_{d_0,k}(\mathbf{p}+\hat{n}s)w_{d_0,k}\right)\mathrm{d}s\\ & =g_y({\mathbf{w}}_{\phi },{\mathbf{w}}_{d_0},\mathbf{\eta }),\end{array}$$

$$\begin{array}{rl}{y}_{t\ast }& =\underset{\mathbf{T}}{\underbrace{\left[\begin{array}{ccc}{n}_{\perp 1}& n_{\perp 2}& 0\\ 0& n_{\perp 1}& n_{\perp 2}\end{array}\right]}}{\mathbf{\varphi }}_{\phi }(\mathbf{x}){\mathbf{w}}_{\phi }=\mathbf{T}({\mathbf{n}}_{\perp }){\mathbf{\varphi }}_{\phi }({\mathbf{x}}_b){\mathbf{w}}_{\phi }\\ & =g_t(\mathbf{w},{\mathbf{x}}_b,{\mathbf{n}}_{\perp }).\end{array}$$

$$\begin{array}{r}p(\mathbf{w})\sim \mathcal{N}(\mathbf{\mu },{\mathbf{\Sigma }}_p)\text{and}p(\mathbf{Y}|\mathbf{w})\sim \mathcal{N}(\mathbf{Y}|g(\mathbf{w}),{\mathbf{\Sigma }}_n),\end{array}$$

$$\log p(\mathbf{Y})\ge \mathbb{E}\left[\log p(\mathbf{Y}|\mathbf{W})\right]-KL\left[q(\mathbf{w})||p(\mathbf{w}|\mathbf{Y})\right]=\mathcal{F}$$

$$\begin{array}{rl}\mathbb{E}\left[\log p(\mathbf{Y}|\mathbf{W})\right]& =\frac{1}{2}\left[N\log 2\pi +\log |{\mathbf{\Sigma }}_n|+(\mathbf{Y}-\mathbb{E}\left[g(\mathbf{w})\right]{)}^{\mathsf{T}}{\mathbf{\Sigma }}_n^{-1}(\mathbf{Y}-\mathbb{E}\left[g(\mathbf{w})\right])\right],\\ KL\left[q(\mathbf{w})||p(\mathbf{w}|\mathbf{Y})\right]& =\frac{1}{2}\left[\text{tr}({\mathbf{\Sigma }}_p^{-1}\mathbf{C})+{\left(\mathbf{\mu }-\hat{\mathbf{w}}\right)}^{\mathsf{T}}{\mathbf{\Sigma }}_p^{-1}\left(\mathbf{\mu }-\hat{\mathbf{w}}\right)-\log |\mathbf{C}|+\log |{\mathbf{\Sigma }}_p|-N\right],\end{array}$$

$$\begin{split}\mathcal{F}&\approx-\frac{1}{2}\big{[}N\log 2\pi+\log|\boldsymbol{\Sigma}_{n}|-\log|\mathbf{C}|+\log|\boldsymbol{\Sigma}_{p}|+(\mathbf{Y}-g(\hat{\mathbf{w}}))^{\mathsf{T}}\boldsymbol{\Sigma}_{n}^{-1}(\mathbf{Y}-g(\hat{\mathbf{w}}))+\left(\boldsymbol{\mu}-\hat{\mathbf{w}}\right)^{\mathsf{T}}\boldsymbol{\Sigma}_{p}^{-1}\left(\boldsymbol{\mu}-\hat{\mathbf{w}}\right)\big{]}\end{split}$$

$$\begin{array}{rl}\mathbf{g}& ={\mathbf{J}}^{\mathsf{T}}{\mathbf{\Sigma }}_n^{-1}(\mathbf{Y}-g(\hat{\mathbf{w}}))-{\mathbf{\Sigma }}_p^{-1}\hat{\mathbf{w}}\\ \mathbf{H}& =-{\mathbf{J}}^{\mathsf{T}}{\mathbf{\Sigma }}_n^{-1}\mathbf{J}+\frac{\partial {\mathbf{J}}^{\mathsf{T}}}{\partial \mathbf{w}}{\Sigma }_n^{-1}(\mathbf{Y}-g(\hat{\mathbf{w}}))-{\mathbf{\Sigma }}_p^{-1}\end{array}$$

$${\hat{\mathbf{w}}}_{k+1}=(1-\alpha ){\hat{\mathbf{w}}}_k+\alpha \mathbf{q}+\alpha \mathbf{\mu },$$

$$\mathbf{C}={\left[{\mathbf{\Sigma }}_p^{-1}+{\mathbf{J}}^{\mathsf{T}}{\Sigma }_n^{-1}\mathbf{J}\right]}^{-1}.$$

$$\begin{array}{rl}\left[\begin{array}{c}{\hat{d}}_0(\mathbf{x})\\ \widehat{\bar{\mathbf{ϵ}}}\end{array}\right]& =\left[\begin{array}{cc}{\mathbf{\varphi }}_{ϵ}(\mathbf{x})& 0\\ 0& {\mathbf{\varphi }}_{ϵ}(\mathbf{x})\end{array}\right]\hat{\mathbf{w}},\\ \hat{\Sigma }& =\left[\begin{array}{cc}{\mathbf{\varphi }}_{ϵ}(\mathbf{x})& 0\\ 0& {\mathbf{\varphi }}_{ϵ}(\mathbf{x})\end{array}\right]\mathbf{C}{\left[\begin{array}{cc}{\mathbf{\varphi }}_{ϵ}(\mathbf{x})& 0\\ 0& {\mathbf{\varphi }}_{ϵ}(\mathbf{x})\end{array}\right]}^{\mathsf{T}},\end{array}$$

$$\mathcal{E}(\mathbf{x})=\left[\begin{array}{c}\frac{P}{EI}(L-x)y\\ -\frac{(1+\nu )P}{2EI}\left({\left(\frac{h}{2}\right)}^2-y^2\right)\\ -\frac{\nu P}{EI}(L-x)y\end{array}\right],$$

$$d_0(\mathbf{x})=c_0\exp \left(-\frac{1}{2}(x-c_1{)}^2/c_2^2-\frac{1}{2}(y-c_3{)}^2/c_4^2\right)+c_5,$$

$${\varphi }_{ϵ,j}=\left[\begin{array}{c}\frac{{\partial }^2}{\partial y^2}-\nu \frac{{\partial }^2}{\partial x^2}\\ -(1+\nu )\frac{{\partial }^2}{\partial x\partial y}\\ \frac{{\partial }^2}{\partial x^2}-\nu \frac{{\partial }^2}{\partial y^2}\end{array}\right]{\varphi }_{\phi ,j}.$$

$$\begin{array}{rl}\frac{{\partial }^2}{\partial x^2}{\varphi }_{\phi ,j}& =\frac{-{\lambda }_{\phi x,j}^2}{\sqrt{L_{\phi x}{L}_{\phi y}}}\sin ({\lambda }_{\phi x,j}(x+L_{\phi x}))\sin ({\lambda }_{\phi y,j}(y+L_{\phi y})),\\ \frac{{\partial }^2}{\partial y^2}{\varphi }_{\phi ,j}& =\frac{-{\lambda }_{\phi y,j}^2}{\sqrt{L_{\phi x}{L}_{\phi y}}}\sin ({\lambda }_{\phi x,j}(x+L_{\phi x}))\sin ({\lambda }_{\phi y,j}(y+L_{\phi y})),\\ \frac{{\partial }^2}{\partial x\partial y}{\varphi }_{\phi ,j}& =\frac{{\lambda }_{\phi x,j}{\lambda }_{\phi y,j}}{\sqrt{L_{\phi x}{L}_{\phi y}}}\cos ({\lambda }_{\phi x,j}(x+L_{\phi x}))\cos ({\lambda }_{\phi y,j}(y+L_{\phi y})).\end{array}$$

$$y_{\ast }=\frac{1}{L}\underset{0}{\overset{L}{\int }}\bar{\mathbf{n}}\left(\sum _j\sum _k{\varphi }_{ϵ,j}(\mathbf{p}+\hat{\mathbf{n}}s)w_{\phi ,j}{\varphi }_{d_0,k}(\mathbf{p}+\hat{n}s)w_{d_0,k}\right)+\left(\sum _k{\varphi }_{d_0,k}(\mathbf{p}+\hat{n}s)w_{d_0,k}\right)\mathrm{d}s$$

$$\begin{array}{rl}& \underset{0}{\overset{L}{\int }}{\varphi }_{d_0,k}(\mathbf{p}+\hat{\mathbf{n}}s)\mathrm{d}s,\\ & \underset{0}{\overset{L}{\int }}{\varphi }_{d_0,k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{{\partial }^2}{\partial x^2}{\varphi }_{\phi ,j}(\mathbf{p}+\hat{\mathbf{n}}s)\right)\mathrm{d}s,\\ & \underset{0}{\overset{L}{\int }}{\varphi }_{d_0,k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{{\partial }^2}{\partial y^2}{\varphi }_{\phi ,j}(\mathbf{p}+\hat{\mathbf{n}}s)\right)\mathrm{d}s,\\ & \underset{0}{\overset{L}{\int }}{\varphi }_{d_0,k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{{\partial }^2}{\partial x\partial x}{\varphi }_{\phi ,j}(\mathbf{p}+\hat{\mathbf{n}}s)\right)\mathrm{d}s,\end{array}$$

$$\begin{array}{rl}{\alpha }_{\phi x}& ={\lambda }_{\phi x,j}(x_0+{\hat{n}}_{1}s+L_{\phi x}),{\alpha }_{\phi y}={\lambda }_{\phi y,j}(y_0+{\hat{n}}_{2}s+L_{\phi y}),\\ {\alpha }_{d_{0}x}& ={\lambda }_{d_{0}x,k}(x_0+{\hat{n}}_{1}s+L_{d_{0}x}),{\alpha }_{d_{0}y}={\lambda }_{d_{0}y,k}(y_0+{\hat{n}}_{2}s+L_{d_{0}y}).\end{array}$$

Full text

1 Introduction

Energy-resolved neutron transmission methods can generate lower dimensional (one- or two-dimensional) images of strain from a higher dimensional (two- or three-dimensional) strain field within a polycrystalline material. The ‘tomographic’ reconstruction of an unknown strain field from these images can be used to study the residual strain and stress within engineering components. Residual stresses are those which remain after applied loads are removed (e.g. due to heat treatment, plastic deformation, etc.), and may have significant and unintended impact on a component’s effective strength and service life — in particular its fatigue life. Measuring and quantifying these strains is important for the validation of predictive design tools, such as Finite Element Analysis, and to aid the development of novel manufacturing techniques — i.e. additive manufacturing.

These strain images are generated by analysing features known as Bragg-edges in the relative transmission of a neutron pulse through a sample. Bragg-edges are sudden increases in the intensity as a function of wavelength and occur when the scattering angle $2\vartheta$ reaches $180^{\circ}$, beyond which no further coherent scattering can occur. The wavelength $\lambda$ at which these Bragg-edges occur can be related to the lattice spacing $d$ within the sample through Bragg’s law: $\lambda=2d\sin\vartheta$. Assuming minimal texture, this can be used to provide a relative measure of strain;

[TABLE]

where $d_{0}$ is the stress-free lattice spacing and $\langle\epsilon\rangle$ is a through thickness average of the normal, elastic strain in the direction of the beam.

The determination of $d_{0}$ is a problem inherent to diffraction and transmission strain analysis. For specific cases where the loading mechanism does not result in changes to the stress-free lattice parameter, its value can be measured prior to loading and in the simplest case (e.g. for an annealed sample) a constant value throughout the sample can be assumed. Several algorithms for strain tomography assuming a known, constant stress-free lattice spacing have been developed. Reconstruction of axisymmetric strain fields is considered in (Abbey et al., 2009, 2012; Kirkwood et al., 2015; Gregg et al., 2017) and more general two-dimensional strain fields in (Gregg et al., 2018; Jidling et al., 2018; Hendriks et al., 2018a).

Many manufacturing techniques (e.g. welding and additive manufacturing) can alter the lattice spacing; for example, as a result of inhomogeneously distributed phase changes (such as the Martensite transformation), or due to gradients in composition as a result of differing chemical states in the starting materials. Since the lattice spacing (in this case $d_{0}$) is sensitive to crystal structure and composition changes, the stress-free lattice parameter may vary throughout the sample. Ignoring variations in $d_{0}$ would cause severe degradation in the quality of a reconstructed strain field. In such cases, measuring $d_{0}$ is more challenging and has been achieved in neutron diffraction measurements by measuring additional directions of strain (Luzin et al., 2011; Choi et al., 2007) and by destructive methods where the strain is relieved by wire cutting the sample into a grid allowing the stress-free lattice spacing to be measured throughout the sample (Paradowska et al., 2005). Although the latter of these two options could be applied to strain tomography it requires the destruction of the sample and creates an additional tomography problem, requiring another set of measurements to be acquired.

Here, we present a method capable of jointly reconstructing the strain field and the $d_{0}$ field from a single set of neutron transmission images. To achieve this both the strain and $d_{0}$ are modelled by a Gaussian process (see for example Rasmussen and Williams (2006)) and equilibrium and boundary conditions are built into the strain model (Jidling et al., 2017). This extends the Gaussian process approach presented by Jidling et al. (2018); Hendriks et al. (2018a) to handle the inherently non-linear nature of this problem. A numerically tractable algorithm based on variational inference (see for example Blei et al. (2017); Jordan et al. (1999)) is provided and the method is validated on a simulated data set.

2 Problem Statement

This paper focuses on the joint reconstruction of the strain field $\boldsymbol{\epsilon}(\mathbf{x})$ and a non-constant stress-free lattice parameter $d_{0}(\mathbf{x})$ from a set of neutron transmission images. Restricting the problem to two dimensions, gives the strain field as the symmetric tensor

[TABLE]

where $\mathbf{x}=\begin{bmatrix}x&y\end{bmatrix}^{\mathsf{T}}$. For brevity, the unique components of strain will be written as $\bar{\epsilon}=\begin{bmatrix}\epsilon_{xx}&\epsilon_{xy}&\epsilon_{yy}\end{bmatrix}^{\mathsf{T}}$ with the coordinate $\mathbf{x}$ omitted where appropriate.

Here, we consider the lattice spacings $d$ as the measurements rather than the standard approach which considers the relative strain of the form form (1). This allows the measurements to be explicitly related to both the strain and the stress-free lattice parameter through the Longitudinal Ray Transform (LRT) (Lionheart and Withers, 2015):

[TABLE]

where $e\sim\mathcal{N}(0,\sigma_{n}^{2})$ and the geometry of each measurement is given by the parameter set $\boldsymbol{\eta}=\{\hat{\mathbf{n}},L,\mathbf{p}\}$; with $\hat{\mathbf{n}}$ as the beam direction, $L$ as the irradiation length, $\mathbf{p}=\begin{bmatrix}x_{0}&y_{0}\end{bmatrix}^{\mathsf{T}}$ as the point of initial intersection between the ray and the sample, and $\bar{\mathbf{n}}=\begin{bmatrix}\hat{n}_{1}^{2}&2\hat{n}_{1}\hat{n}_{2}&\hat{n}_{2}^{2}\end{bmatrix}$. See Figure 1 for the measurement geometry. These measurements are a non-linear function of the two phenomena we wish to estimate; $\boldsymbol{\epsilon}$ and $d_{0}$.

For details on the analysis of neutron transmission data to determine these lattice spacings the reader is referred to Santisteban et al. (2002b, a); Tremsin et al. (2011, 2012). It is also worth noting that the standard deviation $\sigma_{n}$ of these measurements is available.

Furthermore, the strain field inside a sample is a physical property and as such it is subject to equilibrium and boundary conditions (Sadd, 2009). Therefore, it is natural to constrain estimates of the strain field to satisfy these conditions. Using Hooke’s law the equilibrium conditions can be written directly in terms of strain. In two dimensions, this relies on an assumption of plane strain or plane stress. Plane stress is assumed for the remainder of this work, giving the equilibrium conditions as

[TABLE]

where $\nu$ is Poisson’s ratio.

Boundary conditions, in particular the load free surfaces, may also be known. For an unloaded surface, the distribution of forces known as tractions will be zero. Through equilibrium this places additional linear constraints on the strain field, which, assuming plane stress, can be written as

[TABLE]

where $\mathbf{x}_{b}$ is a point on an unloaded surface and $\mathbf{n}_{\perp}$ is the normal to the surface at this point.

An approach to enforcing equilibrium in the estimated strain field is to define a Gaussian process for the Airys stress function from which strain can be derived (Jidling et al., 2018). This non-parametric approach was demonstrated experimentally by Jidling et al. (2018) and compared to other parametric approaches by Hendriks et al. (2018a) with promising results. Boundary conditions in the form of (5) can be included in the estimation process as artificial measurements of zero traction (Hendriks et al., 2018a).

We wish to extend this approach so that both the strain field and the stress-free lattice spacing can be estimated. As the measurements are a non-linear function of the unknowns we cannot directly apply the standard Gaussian process regression methods (Rasmussen and Williams, 2006). There exists several approaches to approximate Gaussian processes for non-linear functions; the Laplace approximation (Rasmussen and Williams, 2006; Bishop et al., 1995), GP variational inference (Steinberg and Bonilla, 2014), and Markov Chain Monte Carlo methods (such as Elliptical Slice Sampling (Murray and Adams, 2010)). For these methods, the measurements are modelled as non-linear functions of the GP sampled at the measurement locations (known as latent function values). The latent function values that best111For a given criterion of best fit, whether it be marginal log likelihood, cross-validation, etc. match the data are determined by one of the above methods. Then, Gaussian process regression is applied with the latent function values taking the place of measurements to determine the function values at the new locations of interest.

The non-static nature of the integral measurement model (3) makes it unclear how to express the measurements as a function of a finite set of latent function values, and hence the above approaches to approximating the GP for non-linear measurements cannot be applied directly. In the following section, we utilise an finite basis function approximation to the GP, and by viewing the problem from an alternate perspective we show how variational inference can be used to solve this non-linear problem.

3 Method

The method presented here is to define a Gaussian process model for the strain field and the stress-free lattice spacing. This Gaussian process model is then approximated using a Hilbert space approximation (Solin and Särkkä, 2014; Jidling et al., 2018). This has two benefits; firstly it removes the need for numerical integration of the covariance function (as discussed by Jidling et al. (2018, 2017)), and secondly it allows us to reformulate the problem as a set of basis functions with unknown coefficients. Variational inference can then be used to learn the coefficients from the LRT measurements and artificial measurements of zero traction.

3.1 Gaussian Process model

The Gaussian process (GP) is a Gaussian distribution of spatially correlated functions;

[TABLE]

The characteristics of the functions belonging to this distribution are governed by a mean function $m(\mathbf{x})$ and a covariance function $k(\mathbf{x},\mathbf{x}^{\prime})$. The covariance function describes the correlation between the function values $f(\mathbf{x})$ and $f(\mathbf{x}^{\prime})$ at any two points $\mathbf{x}$ and $\mathbf{x}^{\prime}$. Careful design of the covariance function can ensure that only functions satisfying desired characteristics belong to the distribution.

Here, we wish to design the covariance function to ensure that only strain fields satisfying equilibrium are contained in the GP prior. Following the formulation in (Jidling et al., 2018) a GP model for the Airys stress functions is defined; $\varphi(\mathbf{x})\sim\mathcal{GP}\left(0,k_{\varphi}(\mathbf{x},\mathbf{x}^{\prime})\right)$. Under the assumption that the sample is plane stress, isotropic, and contiguous, the Airy’s stress functions can be related to strain through the mapping

[TABLE]

where $\mathcal{V}$ is a linear operator, and the superscript denotes which variable the operator acts on. As GPs are closed under linear operators (Papoulis and Pillai, 2002; Hennig and Kiefel, 2013; Wahlström, 2015; Jidling et al., 2018) a GP model for strain that satisfies equilibrium can now be defined;

[TABLE]

where a prior mean function of zero has been chosen.

Additionally, $d_{0}$ function is also modelled by a GP;

[TABLE]

where the prior mean $\bar{d}_{0}$ is chosen to be close to the expected theoretical stress-free lattice spacing for the material used or a measured average in a stress-free sample. The choice of prior mean function does not mean that we believe the $d_{0}$ and $\bar{\boldsymbol{\epsilon}}$ functions to be a particular value, but rather that we do not have any information to suggest otherwise. After the inclusion of measurement information, the mean of the posterior estimate will be updated.

There exists a number of options for the base covariance functions $k_{\varphi}(\mathbf{x},\mathbf{x}^{\prime})$ and $k_{d_{0}}(\mathbf{x},\mathbf{x}^{\prime})$, with both the squared-exponential and the Matérn covariance functions having been successfully used for strain estimation (Jidling et al., 2018; Hendriks et al., 2018a). For a more thorough discourse on available covariance functions the reader is referred to Rasmussen and Williams (2006).

Having defined suitable GP models for the strain and $d_{0}$ fields we now wish to estimate these fields from the LRT and traction measurements. However, the LRT is a non-linear function of these fields and consequently a closed form solution does not exist. The following presents a method for obtaining these estimates that approximates the GP by a finite number of basis functions allowing variational inference to be applied.

3.2 Hilbert Space Approximation to the GP Prior

Here, we make use of the approximation method proposed by (Solin and Särkkä, 2014) and demonstrated to be suitable for the problem of strain tomography (Jidling et al., 2018). This method approximates our covariance function by a finite sum of $m$ basis functions;

[TABLE]

where $S$ is the spectral density of the covariance function. For a stationary covariance function $k=k(\mathbf{r})$, where $\mathbf{r}=\mathbf{x}-\mathbf{x}^{\prime}$, the spectral density and the basis functions are given by;

[TABLE]

where $L_{x}$ and $L_{y}$ control the domain size, and $\boldsymbol{\lambda}=[\lambda_{x},\lambda_{y}]^{\mathsf{T}}$ encodes spatial frequencies of the basis functions. The basis functions are chosen as a solution to the Dirichlet boundary conditions on a rectangular domain, which is a natural choice for the Laplace eigenvalue problem that needs to be solved to approximate the GP (Jidling et al., 2018). The parameters $\theta=\{l_{x},l_{y},\sigma_{f}\}$ are commonly called ‘hyperparameters’ and can be chosen by optimisation (as discussed in Section 5.2). For our application the domain size and spatial frequencies are chosen such that the basis functions spanned a region where their spectral densities, were greater than a minimum threshold. This helps to ensure that the dominant frequencies of the response are captured while maintaining numerical stability.

At this stage, the alternative view point of Bayesian linear regression can be taken. This approach models the unknown function by a set of basis functions with Gaussian coefficients;

[TABLE]

where $\boldsymbol{\phi}(\mathbf{x})$ and $\mathbf{w}$ have dimensions $[1,m]$ and $[m,1]$, respectively. This gives the following model for the strain field $\bar{\boldsymbol{\epsilon}}(\mathbf{x})$ and the stress-free lattice spacing $d_{0}(\mathbf{x})$;

[TABLE]

where the unknown coefficients are independently normally distributed to approximate our GP model; $w_{\varphi,j}\sim\mathcal{N}(0,S_{\varphi}(\boldsymbol{\lambda}_{\varphi,j}))$ and $w_{d_{0},j}\sim\mathcal{N}(\mu_{d_{0},j},S_{d_{0}}(\boldsymbol{\lambda}_{d_{0},j}))$. Where the means $\mu_{d_{0},j}$ are chosen so that the prior has the constant value $\bar{d}_{0}$. In this work, basis functions and parameters corresponding to the $d_{0}$ field will be denoted by the subscript $d_{0}$ and the subscript $k$ will be used as an index. Likewise, basis functions and parameters corresponding to the Airys stress function will be denoted by the subscript $d_{0}$ and the subscript $j$ will be used as an index. The expanded expressions for $\boldsymbol{\phi}_{\epsilon}$ are given in Appendix A.

Using the LRT (3) we can write a model for a predicted measurement as a non-linear function of the unknown coefficients;

[TABLE]

where we have restricted ourselves to a single measurement to simplify the notation. These integrals can be analytically evaluated and the equations are given in Appendix A. Predictions of the boundary tractions $\mathbf{y}_{t}$ at a boundary location $\mathbf{x}_{b}$ with surface normal $\mathbf{n}_{\perp}$ can be written as a linear function of the unknown coefficients;

[TABLE]

The coefficients $\mathbf{w}_{\varphi}$ and $\mathbf{w}_{d_{0}}$ are random variables; as such the predictions $\bar{\boldsymbol{\epsilon}}_{*}$, $y_{*}$, and $y_{t*}$ are also random variables. The problem now is to determine the distribution of the coefficients given a set of LRT and traction measurements. This problem is now in a form allowing variational inference to be used to approximate a solution to the non-linear problem.

3.3 Variational Inference

Variational inference (Blei et al., 2017; Jordan et al., 1999) provides an approximation to the posterior distribution by assuming that it has a certain functional form that contain unknown parameters. These unknown parameters are found using optimization, where some distance measure is minimized. We will in this section provide the details enabling the use of variational inference in solving our problem.

Given $n$ transmission measurements and $n_{t}$ traction measurements, such that a vector of all measurements is given by $\mathbf{Y}=[y_{1},\dots,y_{n},y_{t,1},\dots,y_{t,n_{t}}]^{\mathsf{T}}$, the problem can be written as having prior and likelihood

[TABLE]

where $\mathbf{w}=\begin{bmatrix}\mathbf{w}_{d_{0}}^{\mathsf{T}}&\mathbf{w}_{\varphi}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$, $\boldsymbol{\mu}$ is a vector of all the prior means, and $\boldsymbol{\Sigma}_{p}$ is a matrix with the coefficients prior variance on the diagonals. Here, $g(\cdot)$ is the combined measurement model that expresses the measurement vector $\mathbf{Y}$ as a function of the coefficients. This function is constructed using both (14) and (15). Finally, $\boldsymbol{\Sigma}_{n}=\text{diag}(\sigma_{n}^{2}I_{n\times n},\sigma_{t}^{2}I_{n_{t}\times n_{t}})$, where $\sigma_{t}^{2}$ is a small variance placed on the artificial traction measurements added for numerical reasons.

The non-linear measurement function $g(\cdot)$ makes the likelihood intractable as the prior and likelihood are no longer conjugate. Consequently, the posterior $p(\mathbf{w}|\mathbf{Y})$ is also intractable and so we find an approximate solution using variational inference (Jordan et al., 1999). The idea is to approximate the true posterior by the Gaussian distribution $q(\mathbf{w})\sim\mathcal{N}(\hat{\mathbf{w}},\mathbf{C})$, and find the mean $\hat{\mathbf{w}}$ and covariance $\mathbf{C}$ for this distribution that maximise the Free Energy $\mathcal{F}$. The Free Energy places a lower bound on the log marginal likelihood and hence provides a measure of how well our posterior fits the data;

[TABLE]

where, in this case, $\mathbb{E}[\cdot]$ is the expected value with respect to the approximate posterior $q(\mathbf{w})$ and $KL[\cdot]$ is the Kullback Leibler divergence which provides a measure of difference between the approximate posterior and the true posterior. These terms can be evaluated as (Steinberg and Bonilla, 2014);

[TABLE]

where $N=n+n_{t}$. Here, the expectation of the non-linear function $\mathbb{E}\left[g(\mathbf{w})\right]$ is intractable (Steinberg and Bonilla, 2014) and so the expected maximum is used $\hat{\mathbf{Y}}=g(\hat{\mathbf{w}})$;

[TABLE]

The optimal posterior mean is chosen to maximise the Free Energy. To perform this optimisation a modified Newton’s method is used where the step direction is $\mathbf{q}=-\mathbf{H}^{-1}\mathbf{g}$ and we can calculate the gradient, $\mathbf{g}$, and Hessian, $\mathbf{H}$, of the cost as

[TABLE]

where $\mathbf{J}=\begin{bmatrix}\frac{\partial\hat{\mathbf{Y}}}{\partial\hat{\mathbf{w}}}^{\mathsf{T}}&\frac{\partial\hat{\mathbf{Y}_{t}}}{\partial\hat{\mathbf{w}}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$, and the derivatives and second derivatives are given in Appendix B. At each iteration we update the coefficients according to

[TABLE]

A backwards line search is used to ensure that $\mathcal{F}$ is increased in each iteration. Once the optimal posterior mean is found, the covariance can be found by setting $\frac{\partial\mathcal{F}}{\partial\mathbf{C}}=0$ and linearising about $\hat{\mathbf{w}}$ (Steinberg and Bonilla, 2014), giving;

[TABLE]

Pseudo-code for an algorithm to find approximate distribution of the coefficients $q(\mathbf{w})\sim\mathcal{N}(\hat{\mathbf{w}},\mathbf{C})$ is given in Algorithm 1. Once the coefficients are found, estimates of the strain and $d_{0}$ fields can be estimated. The approximate poster mean and variance for the strain field and stress-free lattice spacing can be computed as

[TABLE]

where $\hat{\Sigma}$ is the joint covariance of the strain and $d_{0}$ estimates. Next, this method is demonstrated on a set of measurements simulated from a theoretical cantilever beam strain field and an artificial $d_{0}$ field.

4 Simulation Results

The method’s ability to jointly reconstruct the strain field and a $d_{0}$ field is demonstrated using simulated measurements. Reconstructions from measurements simulated through two strain fields is shown; the Saint-Venant approximate strain field for a cantilver beam, and a Finite Element Analysis (FEA) strain field from an in-situ loaded C-shape. Additionally, the consequences of ignoring the $d_{0}$ variation on the reconstruction are shown by using the linear measurement model and Gaussian process regression method presented by Jidling et al. (2018) with the addition of traction constraints as shown in Hendriks et al. (2018a). Matlab code to run both examples can be found on Github (Hendriks, 2019).

4.1 Cantilever Beam Example

The method is first demonstrated for the theoretical Saint-Venant cantilever beam as studied. Assuming plane stress, the Saint-Venant approximation to the strain field is (Beer et al., 2010):

[TABLE]

where the geometry is defined in Figure 2. A synthetic stress-free lattice spacing field is defined by

[TABLE]

with the parameters given by $\{c_{0},c_{1},c_{2},c_{3},c_{4},c_{5}\}=\{0.0168,0,7.5\times 10^{-3},7\times 10^{-3},6\times 10^{-3},4.056\}$. The maximum variation $c_{0}$ from a constant base value, $c_{5}$, was chosen to reflect the possible maximum relative variation due to martensitic phase change in $0.8\%$ carbon steel.

Measurements of the form (3) were simulated for $30$ angles evenly spaced between $0^{\circ}$ and $180^{\circ}$, with $100$ measurements per angle, which is on the conservative side based on past experiments (Gregg et al., 2018; Hendriks et al., 2017). The simulated measurements were corrupted with zero-mean noise of standard deviation $\sigma_{n}=c_{5}\times 10^{-4}$ which is equilvalent to $1\text{\times}{10}^{-4}$ standard deviation in strain representing the typical experimental noise (Hendriks et al., 2017; Gregg et al., 2018)

Fifty zero-traction measurements were added along the top and bottom of the cantilever beam for both the presented method and the linear GP regression method. Results are shown in Figure 3. These results show that the presented method successfully reconstructs both the strain field and the $d_{0}$ field with a relative error222$\text{relative error}=\frac{\text{mean}\left(|\text{true}-\text{estimated}|\right)}{\text{max}\left(|\text{true}|\right)}$ of $0.0057$. By contrast ignoring the presence of a $d_{0}$ variation and using a linear GP regression method yields a drastically degraded strain reconstruction with a relative error of $0.3067$.

4.2 In-situ Loaded C-shape Sample

The method is now demonstrated on a more complex strain field given by FEA of a mild steel C-shape sample with geometry defined in Figure 4. The sample was subjected to a $7\text{\,}\mathrm{kN}$ compressive load distributed over $5^{\circ}$ arcs and plane stress was assumed for the analysis. The resulting FEA strain field is shown in Figure 5. This sample and loading conditions correspond to the experimental setup used by Hendriks et al. (2017).

Measurements of the form (3) were simulated through the this strain and a synthetic, smoothly changing $d_{0}$ field is again defined by Equation (25) with parameters given by $\{c_{0},c_{1},c_{2},c_{3},c_{4},c_{5}\}=\{0.01,$5\text{\times}{10}^{-3}$,7.5\times 10^{-3},7\times 10^{-3},6\times 10^{-3},4.056\}$. A total of 60 strain images were simulated with angles evenly spaced between $0^{\circ}$ and $180^{\circ}$, and 180 measurements per image. The simulated measurements were corrupted with zero-mean noise of standard deviation $\sigma_{n}=c_{5}\times 10^{-4}$ which is equivalent to $1\text{\times}{10}^{-4}$ standard deviation in strain representing the typical experimental noise. A total of $131$ zero-traction measurements were added around the boundary of the C-shape excluding the regions within $10^{\circ}$ of the loading points. Reconstruction from the LRT and traction measurements was performed using both the presented method and the linear GP regression method, and the results are shown in Figure 5.

The presented method achieves a mean relative error of $0.023$ and it can be seen that the reconstruction has achieved the correct shape. Whereas assuming a constant $d_{0}$ value gives a mean relative error of $0.137$ and the resulting strain fields show incorrect concentrated peaks in the strain field and areas of tension and compression that are reversed. Despite this improvement there is still some observable difference. In particular, the presented method has a concentrated tensile region on the top left boundary of the C, and does not capture the very concentrated peaks in magnitude on the inside of the C. These peak strains on the boundary are the hardest for the algorithm to reconstruct as they are poorly sampled by the LRT; i.e. they make up only a very small part of each line integral. Additionally, some of this remaining error is due to systematic error in the simulation of the measurements. Which are generated by numerically performing a line integral with each function evaluation being given by an interpolation of the FEA results.

5 Additional Remarks

5.1 Sensitivity to the Traction Measurement Variance

Boundary conditions given by unloaded surfaces are a natural inclusion as they are an artefact of the physical world. This information is included in the form of artificial measurements of zero traction, however it was found that a small variance needed to be placed on these measurements and the rate of convergence was impacted by the size of this variance. Conceptually, this variance is analogous to a constraint tolerance for optimisation procedures. Too large a variance (or not enough traction measurements) and the algorithm may fail to converge to the correct strain field. This indicates that the traction measurements are ensuring that the problem is observable, which is supported by the findings of Hendriks et al. (2018b) where the inclusion of traction measurements allowed a constant $d_{0}$ value to be found as a hyperparameter. Conversely, too small a variance and the algorithm is unable to take optimisation steps of significant size, resulting in a large number of iterations to converge. Methods for optimally choosing this variance is an avenue for future research.

Despite this, it was found that the algorithm worked well over a reasonable range of traction variances. Typically the standard deviation of the traction measurements could be set two orders of magnitude smaller than the measurement standard deviation or in the range of $1\times 10^{-5}$ to $1\times 10^{-7}$.

5.2 Hyperparameter Optimisation

The hyperparameters $\theta=\{l_{x},l_{y},\sigma_{f}\}$ can be found by performing an optimisation using $\mathcal{F}$ as the objective function. However, the gradients of $\mathcal{F}$ with respect to the GP hyperparameters are not trivial and so Steinberg and Bonilla (2014) suggests that gradient free optimisation approaches could be used. In this work, both Bayesian optimisation (Mockus, 2012) and the Nelda-Mead method (Nelder and Mead, 1965) were found to work; with the Nelda-Mead method requiring less computation time.

5.3 Computational cost

The majority of the computational burden comes from building the building the forward model of the LRT measurements. In particular, all the combinations of the basis function of the stress-free lattice spacing and the strain field need to be computed for each measurement. The resulting matrix has $nm_{d_{0}}m_{\varphi}$ elements. However, as this matrix only needs to be computed once, it is still feasible to solve even with a large number of basis functions.

6 Conclusion

This paper considers an extension of the strain tomography problem where the stress-free lattice parameter is a known constant, to the more general case where it is unknown and varies throughout the sample. A method for the joint reconstruction of a strain field and a varying stress-free lattice parameter from a set of neutron transmission strain images has been presented. This method extends the Gaussian process based approach previously used for strain tomography to the subsequently non-linear problem, and ensures that the estimated strain fields satisfy equilibrium and can include knowledge of boundary conditions. This was achieved by reformulating the problem in terms of basis functions and unknown coefficients. Variational inference was then employed to find estimates of the coefficients.

The method was tested on a set of simulated data, and importantly, these results demonstrate that it is possible to perform this joint reconstruction. Further, the results obtained by ignoring variations in $d_{0}$ and applying the linear GP regression method are provided and show that this assumption, if incorrect, severely degrades the accuracy of reconstruction.

Future work will involve planning an experiment to acquire a data set on which to further evaluate the methods performance.

Acknowledgements

This work is supported by the Australian Research Council through the Discovery Project scheme (ARCDP170102324), as well as the Swedish Foundation for Strategic Research (SSF) via the project ASSEMBLE (contract number: RIT15-0012), and by the Swedish Research Council via the projects Learning flexible models for nonlinear dynamics (contract number: 2017-03807) and NewLEADS - New Directions in Learning Dynamical Systems (contract number: 621-2016-06079).

Appendix A Basis Functions

The basis functions for the strain field were defined in (13) as

[TABLE]

As such, the components of $\phi_{\epsilon,j}$ can be built from

[TABLE]

A predicted LRT measurement was defined by (14) as

[TABLE]

where for clarity we restrict ourselves to a single measurement. Therefore, we need the components

[TABLE]

To make the expressions briefer, we introduce the notation

[TABLE]

Giving

[TABLE]

Appendix B Measurement Model Derivatives

Here we give the derivates of the measurement model $g(\cdot)$ about the current $\hat{\mathbf{w}}$. The measurement model is a concatenation of equations (14) and (15), and so we require the derivatives of the predicted LRT measurement, $\hat{y}$, and the predicted traction $\hat{y}_{t}$. For clarity we restrict the following to a single $\hat{y}$ and $\hat{y}_{t}$. The first derivatives are given by

[TABLE]

The second derivates are

[TABLE]

where

[TABLE]

Explicit formulation of the second derivatives allows the cost functions curvature to be taken into account in the optimisation procedure, greatly improving the rate of convergence.