Energy Resolved Neutron Imaging for Strain Reconstruction Using the Finite Element Method

A novel pulsed neutron imaging technique based on the finite element method is used to reconstruct the residual strain within a polycrystalline material from Bragg edge strain images. This technique offers the possibility of a nondestructive analysis of strain fields with a high spatial resolution. The finite element approach used to reconstruct the strain uses the least square method constrained by the conditions of equilibrium. This inclusion of equilibrium makes the problem well-posed. The procedure is developed and verified by validating for a cantilevered beam problem. It is subsequently demonstrated by reconstructing the strain from experimental data for a ring-and-plug sample, measured at the spallation neutron source RADEN at J-PARC in Japan. The reconstruction is validated by comparison with conventional constant wavelength strain measurements on the KOWARI diffractometer at ANSTO in Australia. It is also shown that the addition of a Tikhonov regularisation scheme further improves the reconstruction.

In these experiments, the term Bragg edge [6,7] refers to a sudden increase in the relative transmission of a neutron beam passing through polycrystalline solids as a function of wavelength. A neutron, of wavelength λ, can be coherently scattered by crystal planes with lattice spacing d, provided that the scattering angle θ satisfies Bragg's law (λ = 2d sin θ). A sudden increase in transmission occurs once λ = 2d is exceeded as a neutron cannot be scattered by more than 180 • [3,8], so neutrons are backscattered, and no further diffraction occurs from that particular plane [9,10].
While other approaches exist [2,[11][12][13], the process of measuring Bragg edges we use here relies on the measurement of the transmission spectra using the time-of-flight or energy-resolved techniques [14]. This method requires a pulsed neutron source. Such a pulsed neutron source can be found in Japan (J-PARC) [15,16], UK (ISIS) [17], and USA (SNS). The greatest advantage of neutron strain tomography is that the incident beam flux is fully utilised, helping to reduce the data collection time. Modern technology uses a pixelated detector consisting of an array of up to 512 × 512 pixels with spatial resolution as small as 55 µm [18]. Such strain imaging raises the prospect of strain tomography, and several attempts have been made to solve the resulting tensor reconstruction problem over the past decade. Often the methods of reconstruction have revolved around special cases, for example, axial-symmetry [19,20].
It is essential to note that the inverse problem is ill-posed, and reconstruction of the strain is not easy without imposing further conditions [21] or by making simplifying assumptions [22]. The question of how best to impose these extra conditions has been the subject of intense debate, and different approaches have been proposed in the literature. Recently, it has been shown that by applying the condition of either equilibrium or compatibility, the reconstruction is possible [22,23] (we note that the main difference between these works is the choice of basis functions which represent the strain field).
In this paper, we describe a method by which it is possible to tomographically reconstruct the elastic strain from a series of Bragg edge strain measurements using a finite element discretisation constrained by equilibrium. This condition fits naturally into the finite element framework. This method offers several advantages over previous methods due to the very desirable properties of the finite element method. In particular its ability to accurately represent highly spatially varying functions. The proposed algorithm is tested on a cantilevered beam simulated data in two dimensions. It is shown to be capable of reconstructing a strain tensor field after imposing the equilibrium conditions [24,25]. The algorithm is then applied to experimental data for a ring-and-plug geometry. We introduce a smoothing function to the minimisation problem with a regularisation parameter. Hence, minimising the value of the objective function will give us a regularised resistivity update equation to reduce the noise in the reconstructed images.

Longitudinal Ray Transform
We outline here the experimental technique which has recently been developed that provides information on the average strain component in the direction of the incident beam [2]. As mentioned earlier, Bragg edges are formed by backscattering radiation. Hence, relative shifts in their position provide a measure of the average normal strain within a sample in the direction of the beam [26], i.e., Therefore, the average strain within a body as measured by Bragg edge neutron transmission can be idealised as a line integral known as Longitudinal ray transform (LRT) which captures the average component of strain along the line s in the direction of the unit normaln = (n i , n j ) T = (cos ϑ, sin ϑ) T . We define , y(s, a)) n j ds, where ε ij is the component of tensor strain field ε ∈ R 2×2 , which is mapped to an average normal component of a strain in the direction ofn. The ray enters the sample at the position x a and L is the ray length inside the sample for a particular angle ϑ. This configuration is shown in Figure 1. This technique relies on the overall change in the lattice spacing along the ray [21,23]. Measurements are taken in each orientation, ϑ i , where a profile is measured of the form Γ ε (x, y, ϑ i ). While the inherent symmetry of the transform implies 180 • are sufficient, in practice, measurements need to be taken over an entire revolution, i.e., 360 • . Lionheart and Withers [21] demonstrated that the integral line LRT is a non-injective map from ε → Γ ε (x, y, ϑ) and hence the strain field produced by any given set of projection is not unique [27]. As a consequence, it is not possible to reconstruct the strain distribution within a body in the general setting. Hence, additional information (equilibrium or compatibility constraints) is required to ensure it is the actual strain field which is recovered from all the possibilities. To this end, several prior approaches have been developed that rely upon assumptions of compatibility or equilibrium to constrain the problem further. Compatible strain fields are those that can be written as the gradient of a displacement field in a simply connected body (i.e., conservative strain fields.). For Example, Abbey et al. [23] developed an algorithm using different radial basis functions (since their problem was axisymmetric), and they enforced the compatibility constraints [22]. Special cases have been considered including axis-symmetric systems [19,20,23,28] and granular systems [29] with equilibrium conditions imposed. The unknown strain can also be reconstructed by using a machine learning technique known as Gaussian Process [25,30], where equilibrium is used as a central technique to ensure that the strain is chosen uniquely. Furthermore, arbitrary strain fields produced due to in-situ loadings have been reconstructed by using compatibility [22,31]. We present here a reconstruction using the finite element method, noting that this method has proved itself to be widely applicable.

Solution using Finite Element Basis Functions
In our numerical implementation, each component of the strain is approximated by a linear combination of basis functions. These basis functions come from the finite element method. The line integral (2) is solved in terms of the unknown strains, and it is equated to the Bragg edge measurement. The uniqueness of the solution is guaranteed by the equilibrium equation, which is imposed on the minimisation problem used to calculate the backward map for the strain, in the form of extra constraints. We formulate the problem as follows is the symmetric strain tensor field, i.e., ε 12 (x, y) = ε 21 (x, y), is normal component, and L = length of a ray inside the sample/geometry. The main problem is to find the ray transform of the components of strain aligned with the direction of projectionn defined in Equation (2) rewritten in the form where x a = (x a , y a ) and x b = (x b , y b ) are the entry and exit points of the ray respectively. The computational solution of the integral Equation (3) requires discretisation, i.e., the integral is expressed in terms of finitely many unknowns. We discretise the sample by using a quadrilateral mesh with m nodes and P elements. A visualisation of such a field over a rectangular sample discretised into rectangles is shown in the Figure 2. A given ray can enter and exit the sample at arbitrary points. Hence, by applying discretisation we obtain the approximate problem as follows where P = {P 1 , P 2 , . . . , P n } is the set of elements and L P is the length of the ray in each element. Note that this length will be zero in many of the elements. An example of such a discretisation is shown in Figure 2 where the first ray is intersecting with the elements 1, 2 and 5, whereas the second ray is intersecting the elements 2, 3, 5 and 6. The strain in any element depends on the strain value at the corner of the quadrilaterals which are the unknowns. In general each node strain will influence the strain in four elements. The strain at each node has three components. The strain is expressed using the standard basis function as follows, see Equation (5) and (7): where, β P ij , γ P ij , η P ij and ζ P ij are the coefficients which are determined from the nodal values of strain in the element P. The measurement for the first ray in Figure 2 can be approximated by where L = ∑ P L P , L P is the segment of the ray inside the P th rectangle where P = 1, 2 and 5 and ε P ij is strain value in the element P. The value of β P ij , γ P ij , η P ij and ζ P ij for each element can be found by solving the following equation  where x i and y i are the coordinates of the corner points of the quadrilaterals. In general, each element P will have different ray entry and exit points. Using the above line integral expression and evaluating the basis functions for each element, the integral can be reformulated in terms of nodal strain and we obtain a system of equations of the form where N is the total number of projections, Γ v , v = {1, ..., N}, is the value from each measurement, K is the matrix derived from the integrals in Equation (4) expressed in terms of the nodal strain through Equation (7). We can write this in compact form as where Γ is a vector containing all of the Bragg edge strain measurements, K is the coefficient matrix with elements that contain unit direction vector components and shape function evaluations which will be a sparse matrix, and = ε 1 11 · · · ε m

22
T is a vector containing all the unknowns for each element. Once the matrix K and vector Γ have been formed, the problem is reduced to one of solving the linear algebraic system of equations for the unknown coefficients represented by vector . In practice, the system is usually over-determined since the number of unknown coefficients is relatively small compared to the amount of experimental data available. Furthermore, the K matrix will be sparse.
As it was pointed out before by Lionheart and Withers [21] the strain field is not uniquely defined within an object from these measurements. For this reason with apply the constraints to our problem obtained by solving the equilibrium equations. From Hooke's law, the equilibrium equation can directly be written in terms of strain. Assuming plane stress condition, the equilibrium conditions can be written as where, ν is the Poisson's ratio. To reconstruct the strain, the equilibrium Equations (10a) and (10b) are integrated over each element P, which will lead us to the following equations using (5) This provides another set of a system of equations where C represents the equilibrium integral matrix, which has two rows. Solutions to the minimisation problem were found by least-square fitting [32], where the problem is reduced to: find a vector such that min The minimisation problem (13) is solved straightforwardly using least squares.

Cantilevered Beam
To demonstrate the performance of the proposed algorithm, a well-known 2D cantilevered beam problem is studied, which was previously examined by Wensrich et al. [31]. We consider the 2D strain field for beam geometry of the rectangle [0, 12] × [0, 10] with the load P of 2 kN displayed in Figure 3. Material properties of the beam are representative of common steel, whereas other parameters are mentioned below. This beam problem is excellent for testing the algorithm since the analytical solutions to the strain field exist.
Assuming plane stress, the Saint-Venant approximation to the strain field is [33]: where I is the second moment of area, L is the total length of the beam, W is the width of the beam, t is the thickness, E is Young's modulus, ν is Poisson's ratio, and the dimensions are shown in Figure 3.  We found that the proposed reconstruction algorithm is extremely effective in achieving strain field reconstruction. A finite element model of the system was constructed, with a structured quadrilateral mesh size 4 × 4. Simulation results suggest that the reconstruction algorithm can converge to an adequate reconstruction provided that measurements are taken over the entire 360 • of a sample. Problem discretisation and numerical errors can undoubtedly contribute to an imperfect reconstruction (with more noise). Rapid convergence to the true solution was observed as the number of projections was increased. This convergence provides confidence in the ability of the algorithm to converge to a true solution in the presence of real experimental uncertainties.

Reconstruction of the Offset Ring-And-Plug
We now test the algorithm on experimental data for an offset ring-and-plug sample, which was used previously [24]. The sample geometry of the offset ring-and-plug is shown in Figure 5. The described sample contained a total interference of 40 ± 2 µm produced through cylindrical grinding. More details about the sample can be found in [24]. A steel bar EN26 was heated to relieve stress and provide a uniform structure prior to the assembly. The final hardness of the sample was 290 HV. The strain profile was measured on RADEN together with an MCP detector at a distance of 17.9 m from the source of the beam. RADEN is an energy-resolved neutron imaging instrument at the Japan Proton Accelerator Research Complex (J-PARC), Japan, [15,16] was used to obtain the relative shifts of the Bragg edge corresponding to the lattice plane of the offset ring-and-plug steel sample. Neutron strain scanning was carried out on KOWARI, a residual stress diffractometer at the Australian Nuclear Society and Technology Organisation (ANSTO), Australia to provide independent validation of our reconstruction. Sampling times on KOWARI were based on providing uncertainty in strain around 7 × 10 −5 , which required around 30 h of beamtime per component. However, sampling times on RADEN were based on statistical uncertainty of the order 1 × 10 −4 . In total, 50 profiles were measured at golden angle increments in ϑ with a sampling time of 2 h per projection. A finite element quadrilateral mesh is used to discretise the domain of the sample, as shown in Figure 6. Two types of mesh patches have been considered: structured and unstructured meshes. Reconstruction results can be seen with each mesh type in Figures 7 and 8. Again, as mentioned before for the beam problem that the reconstructed strain field contains noise, which can be from numerical discretisation or measurements. Some techniques are explained in the next section to cope up with noise present in the strain field.   Different type of the mesh has been shown variation results, which proves that our algorithm is highly dependent on the mesh. It was observed that results with an unstructured mesh show better agreement than the structured mesh in terms of noise. This is because the unstructured mesh is evenly distributed throughout the sample domain, unlike the structured mesh, resulting in reduced numerical noise. Results are then compared with the pointwise measurement of strain on KOWARI, the constant wavelength diffractometer shown in Figure 9 with the reconstructed transmitted measurements of strain on RADEN shown in Figures 7 and 8.

Tikhonov Regularisation
Until this point our reconstruction does not involve any smoothing. As a result, reconstruction images have noise which can arise from different sources such as from instrumental measurement noise or the simulation procedure. To manage this noise, Tikhonov regularisation is used [34]. The Tikhonov regularised estimate is defined as the solution of the following minimisation problem where the first term is the same euclidean norm used before in Equation (13). The second term is known as the regulariser which captures the prior knowledge and behavior of through an additional penalty term where · 2 is Euclidean norm, α > 0 is the regularisation parameter which specifies the amount of regularisation. The effect of regularisation can be varied due to the scale of the matrix B. The matrix B is a block diagonal matrix where the block diagonal entries can be chosen in several ways such as zero matrix which will bring our problem back to unregularised least square problem, it can be identity matrix shown in Figure 10. Hence, in our case, B is chosen as the block diagonal matrix as where S ij = Ω φ i · φ j for i, j = 1, 2, ..., n, and φ i , ∀ i, are the standard basis functions for quadrilateral. Numerically, the minimum is achieved by solving a linear least-square problem of the form: Above equation is solved in MATLAB with a built-in function "lsqlin". The main problem here is to determine proper regularisation parameter α; if the parameter is significant, the solution will deviate from the correct solution, and if the parameter is small, then there will not be any significant difference in the noise. For now, we find this parameter by trial and error; however, it can also be achieved using optimisation methods [35]. The effect of Tikhonov regularisation can be seen in the Figures 11-13, where the difference is shown for different values of α. Figure 10. Regularised strain field ε xx , ε xy , ε yy respectively for unstructured mesh type, with S as identity matrix and α = 0.005. Figure 11. Regularised strain field ε xx , ε xy , ε yy respectively for unstructured mesh type, with S as stiffness matrix and α = 0.001.

Figure 12.
Regularised strain field ε xx , ε xy , ε yy respectively for unstructured mesh type with S as stiffness matrix and α = 0.005. Figure 13. Regularised strain field ε xx , ε xy , ε yy respectively for structured mesh type, with S as stiffness matrix and α = 0.005.
To summarise, in Tikhonov regularisation, we approximate the minimum norm by least squares. Least square solution φ R depends on Kφ, by a vector depending on the regularisation parameter α > 0. Reconstruction is done on finite element rectangular mesh with 3688 elements for the unstructured mesh and 3776 elements for structured mesh type. Minimisation problem for both cases (with and without regularisation) is solved in MATLAB by using built-in function 'lsqlin', which were then plotted by using 'scatteredInterpolant' with linear map fitting. The proposed algorithm does not solely depend on the sample geometry and hence, can be extended to three-dimensional sample bodies. The true difficulty will be the computational cost and time since the size of the problem will be larger.

Conclusion
The proposed method helps to reconstruct the entire strain field, satisfying equilibrium with no assumptions of compatibility and hence is suitable for reconstructing residual strain fields. The KOWARI strain and RADEN reconstructed strain field measurements show closer agreement with the inclusion of regularisation. The proposed method was validated with simulated data and strain estimates from experimentally measured data at J-PARC, Japan, which were compared to the strain calculation from a conventional diffraction method obtained at KOWARI. In two dimensions, full strain fields tomography using Bragg edge images can now be achieved using physical constraints as equilibrium. This method opens up further research for future investigations, including extending this technique to three dimensions. Furthermore, the proposed method allows us to use adaptive meshes that can focus on the highly strained area in the sample, which can be achieved by calculating a gradient over a sample.